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Abstract 


The  development  and  use  of  Micro  Air  Vehicles  is  becoming  more  and  more 
important  to  the  military.  Size,  weight,  maneuverability,  stealth,  and  fuel  consumption  are 
just  some  of  the  constraints  on  a  future  Micro  Air  Vehicle.  Perfect  blueprints  for  a  small 
scale  flying  machine  are  both  insects  and  humming  birds.  Medium  sized  hummingbirds 
are  found  to  have  wings-beats  at  a  frequency  of  18-28  Hz.  There  is  a  vast  amount  of 
complexity  to  just  how  these  creatures  can  create  lift;  however,  this  study  looks  at 
different  beams  made  of  common  materials  with  1st  natural  bending  frequencies  in  this 
range.  This  study  documents  and  evaluates  the  use  of  analytical  tools  to  solve  for 
nonlinear  characteristics  of  a  system.  A  function  called  the  backbone  curve  is 
incorporated  into  a  MATLAB  program.  The  Hilbert  transform  characterized  by  the 
nonlinear  decrement  approach  captures  all  of  the  necessary  coefficients  for  this  function. 
ABAQUS/CAE  is  relied  upon.  The  numerical  finite  element  results  are  compared  to 
experiments.  This  work  gives  a  better  understanding  of  how  materials  and  geometry 
perform  when  used  in  Micro  Air  Vehicle  design  considering  large  displacements. 


FINITE  ELEMENT  SOLUTION:  NONLINEAR  FLAPPING  BEAMS  FOR  USE 


WITH  MICRO  AIR  VEHICLE  DESIGN 


I.  Introduction  /  Literature  Review 


Chapter  Overview 

The  purpose  of  this  chapter  is  to  introduce  the  need  for  a  Micro  Air  Vehicle,  its 
possible  solutions,  similar  projects  underway,  and  past  relevant  research. 

Background 

Micro  Air  Vehicle  (MAV)  refers  to  a  relatively  new  unmanned  type  of  aircraft 
that  are  generally  an  order  of  magnitude  smaller  than  typical  aircraft  of  today.  A  Micro 
Air  Vehicle  was  defined  by  the  Defense  Advanced  Research  Project  Agency  (DARPA) 
MAV  program  in  1997  as  “vehicles  no  larger  than  six  inches  in  either  length,  width  or 
height  and  perform  a  useful  military  mission  at  an  affordable  cost.”  [17]  These  MAV(s) 
are  foreseen  to  be  expendable  and  available  to  individual  soldiers  or  units  for 
requirements  in  reconnaissance,  surveillance,  battle  damage  assessment,  targeting,  sensor 
placement,  communication  relays,  or  for  sensing  chemical,  nuclear  or  biological 
substances  and  be  able  to  sustain  missions  of  up  to  2  hours  long  [17].  There  are  three 
main  types  of  Micro  Air  Vehicles  that  could  offer  a  solution  to  these  needs:  fixed  wing 
vehicles,  rotorcraft  vehicles  and  flapping  wing  vehicles  or  ‘omithopters’.  Though  all  are 


1 


feasible  solutions,  a  flapping  wing  design  is  of  interest  because  it  allows  best  for  loitering 
in  small  areas  while  also  having  a  small  acoustic  signature. 

The  most  useful  flapping  wing  Micro  Air  Vehicle  design  would  be  one  that 
resembles  the  flight  of  a  hummingbird  for  its  sustained  hovering  capabilities  along  with 
vertical  and  horizontal  movements  that  are  produced  by  the  flapping  of  its  wings.  This 
same  flight  pattern  can  be  found  in  certain  bats  and  insects  as  well.  This  would  allow  a 
MAV  access  to  areas  not  reachable  at  present  times.  The  Giant  Hummingbird  with  a 
length  of  8.5  inches  and  weight  of  20  grams  closely  resembles  the  size  of  an  MAV  and 
has  a  wing  frequency  of  approximately  18-30  Hz  [15].  Smaller  hummingbirds  have  wing 
frequencies  of  60-80  Hz  while  insect  wing  frequencies  are  even  higher  such  as  a  bee 
wing  of  200  Hz.  A  MAV  or  omithopter  with  similar  flapping  frequencies  is  eventually 
desired. 

Omithopter: 

The  term  omithopter  comes  from  the  Greek  word  ornithos  "bird"  and  pteron 
"wing".  Omithopters  are  defined  as  mechanical  aircraft  which  mimic  the  flight  of  birds, 
bats  or  insects.  The  flapping  wing  motion  is  used  for  most  flight  characteristics  including 
lift,  propulsion  and  control.  While  mimicking  the  flapping  wing  motion  of  birds,  bats  or 
insects,  omithopters  are  usually  built  to  the  same  scale  as  these  creatures  as  well,  but 
don’t  have  to  be  [19].  Inspired  by  Leonardo  Da  Vinci,  many  successful  attempts  at 
manned  flight  in  omithopters  date  back  to  1870.  With  modem  technological  advances, 
smaller  scaled  omithopters  are  being  developed  all  over  the  world. 


2 


Relevant  Research 


‘Fly  Tech’  recently  developed  a  widely  commercialized  remote  controlled 
dragonfly  with  two  sets  of  18”  wings  as  seen  in  Figure  1  [20]. 


Figure  1.  Fly  Tech  Remote  Control  Dragonfly 


This  vehicle  shows  that  a  battery  powered  remote  control  omithopter  can  be 
produced  at  an  affordable  cost;  however,  dragonfly  is  three  times  as  large  as  a  defined 
Micro  Air  Vehicle.  There  is  also  no  advanced  control  with  this  model.  The  wings  can 
flap  faster  to  produce  more  lift  and  propulsion  while  turning  is  accomplished  using  a 
small  tail  rotor.  Hovering,  vertical  and  backwards  flight  movements  are  not  replicable  on 
a  regular  basis.  The  battery  has  only  a  10  minute  lifespan  before  recharge  is  necessary. 

In  July  2006  the  University  of  Arizona  placed  first  in  the  Omithopter  portion  of 
the  10th  Annual  Micro  Air  Vehicle  Competition  in  Provo,  Utah  with  their  vehicle  shown 
in  Figure  2  [14].  The  University  of  Florida  continues  to  dominate  in  the  other  three 
portions  of  the  MAV  competition  [1], 
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Figure  2.  University  of  Arizona’s  Ornithopter 


The  ornithopter  competition  is  judged  on  both  size  and  laps  the  aircraft  can 
complete  in  2  minutes.  The  University  of  Arizona’s  omithopter’s  wingspan  measures 
less  than  6  inches  and  is  the  world’s  smallest  radio  controlled  ornithopter  [14].  The 
wings  are  made  of  carbon  fiber  struts  covered  with  Mylar  and  flap  similar  to  a 
hummingbird  at  a  frequency  of  50  Hz.  This  aircraft  is  one  of  the  closest  developments  to 
date  of  what  a  future  flapping  wing  MAV  will  be;  however,  this  ornithopter  only  has  a 
flight  time  of  3  minutes  and  while  it  flies  like  a  bird  it  too  does  not  have  the  hovering, 
vertical,  and  backward  flight  characteristics  of  a  hummingbird. 

A  Micro  Air  Vehicle  and  wing  capable  of  replicating  the  intricacies  of  a 
hummingbird  in  flight  is  still  years  away.  There  are  many  different  size,  shape,  material, 
and  design  possibilities  to  list.  A  fundamental  element  of  an  engineering  structure  is  a 
beam  and  to  simplify  things,  most  structures  can  be  modeled  as  such,  including  bridges, 
antenna  towers,  helicopter  rotors,  and  wings.  Even  a  complex  wing  structure  is  made  up 
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of  multiple  beam  elements  and  the  understanding  of  a  beam  structure  will  benefit  the 
understanding  of  the  more  complex  structure  [13].  Pramod  Malatkar  understood  this  and 
in  July  of  2003  put  in  print  an  extensive  study  on  nonlinear  vibrations  in  metallic 
cantilever  beams  and  plates  subject  to  transverse  harmonic  excitations.  This  study 
includes  parametric  system  identification,  and  jump  frequency  determination,  but 
concentrated  on  energy  transfer  between  widely  spaced  modes.  He  states  that  interesting 
phenomenon  occurs  in  structures  that  exhibit  nonlinearity  that  cannot  be  explained  by 
linear  models  [13].  Most  structures  experience  nonlinear  characteristics  although  some 
are  weaker  than  others.  Linear  models  are  only  accurate  when  displacements  are  kept 
small  or  when  the  nonlinearity  is  weak. 

In  2003,  a  combined  effort  of  Robert  Gordon,  Joseph  Hollkamp,  and  Michael 
Spottswood  did  a  study  on  sonic  fatigue  problems  in  military  aircraft.  “Sonic  fatigue  is 
characterized  by  large  amplitude,  nonlinear  response  of  thin  skins  to  intense  aeroacoustic 
loads”  [9],  The  main  purpose  of  their  research  was  to  set  up  well  characterized 
experiments  that  could  be  replicated  with  Finite  Element  methods  for  validation  of  the 
FEMs  themselves.  These  experiments  were  for  a  9.00”  clamped-clamped  steel  beam. 
One  experiment  recorded  simple  harmonic  response  data  for  the  beam  to  acquire  linear 
frequencies  (coo),  and  the  cubic  stiffness  coefficient  (P)  for  the  beam  using  both  the 
Hilbert  transform  and  the  Estimated  Backbone  equation  [9], 

Also  in  2003,  Matthew  Ernst  and  Anthony  Palazotto  conducted  a  study  of  the 
same  clamped-clamped  Steel  beam  using  ABAQUS/CAE  for  the  simulation  of  their 
simple  harmonic  response  data  [6],  They  applied  both  the  Hilbert  transform  and  the 
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estimated  Backbone  Curve  equation  to  their  results  and  found  the  natural  frequencies  and 
the  values  of  linear  coo  and  the  cubic  stiffness  coefficient  (P)  for  various  values  of  the 
estimated  damping  ratio  (£).  The  results  were  compared  to  those  in  [9]  for  accuracy. 

Yamaki  and  Mora  conducted  an  experimental  investigation  in  1980  of  a  Clamped- 
Clamped  beam  with  harmonic  base  excitation  on  an  electrodynamic  shaker.  FRF  curves 
were  generated  with  swept  sine  excitation  for  several  input  forces  [16]. 

Capt  Adam  Tobias  conducted  a  study  of  both  Aluminum  and  Composite  clamped- 
free  beams  from  2006-2007,  concentrating  on  visual  methods  to  capture  the  beams 
vibrations.  These  methods  include  high  speed  photography,  and  multiple  laser 
vibrometers.  The  beams  in  the  study  included  a  clamped-free  10.125”  Aluminum  beam 
and  a  clamped-free  6.08”  Composite  beam  [15]. 

Summary 

This  study  will  look  at  a  finite  element  analysis  using  ABAQUS  for  five  beams, 
its  ability  to  show  and  provide  data  for  the  nonlinear  characteristics  of  the  beams. 
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II.  Theory 


Linear  Systems 

One  of  the  most  widely  used  simplifications  in  beam  theory  is  the  assumption  that 
a  system  is  linear.  Linear  systems  are  systems  that  are  subject  to  superposition  where 
their  input  and  outputs  can  be  added  together  yielding  another  set  of  correct  results  for 
that  system.  Likewise,  the  inputs  and  outputs  can  be  multiplied  by  constants  as  well. 

This  assumption  not  only  makes  calculations  much  easier,  but  it  is  also  very  accurate  for 
small  deflections.  The  strain  for  an  Euler-Bemoulli  beam  is  shown  by  the  following 
equation  [3]. 


£x  £;txial  T  ^bending 

Saxial  =  du/dx 

lending  =  Z  (d2w/dx2) 

where 

u  =  displacement  in  the  axial  direction 
w  =  displacement  in  the  direction  of  the  force  P 
z  =  distance  from  the  midsurface  of  the  beam 
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For  a  beam  under  loading  conditions  shown  above  in  Figure  3,  the  axial  strain 
term  du/dx  in  linear  theory  will  be  equal  to  zero  as  the  rotation  of  the  beam  is  not  great 
enough  to  add  a  component  of  force  in  the  axial  direction  on  the  beam.  For  a  beam  of 
small  deflections,  these  equations  work  well  in  approximating  the  strain,  however,  in 
everyday  life,  no  physical  system  is  strictly  linear  [13].  The  flapping  motions  of  the 
beams  in  this  study  well  surpass  the  small  deflections  for  a  linear  system.  With 
deflections  much  greater  than  the  thickness  of  the  beam,  the  beams  will  be  better  suited 
for  nonlinear  analysis. 

Nonlinear  Systems 

The  term  nonlinearity  means  that  the  inputs  and  outputs  of  a  system  are  not 
additive  as  they  are  in  a  linear  system.  The  response  is  not  directly  proportional  to  the 
inputs  that  produced  it  [13].  There  are  three  main  types  of  nonlinearities  in  a  system, 
material  nonlinearity,  contact  nonlinearity,  and  geometric  nonlinearity.  One  example  of 
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geometric  nonlinearity  stems  from  deformations  large  enough  where  the  system  develops 
membrane  stretching  forces  that  carry  a  portion  of  the  load  [3].  These  new  forces  bring 
about  membrane  strain  terms  that  must  be  accounted  for  in  our  nonlinear  beam  strain 
equations.  Large  rotations  are  also  considered  and  in  turn,  the  strain  due  to  beam  bending 
will  change  as  well. 

£x  —  ^membrane  T  ^bending  (2) 

2  2 

^membrane  =  du/dx  +  V2  (du/dx)  +  V2  (dw/dx) 

^bending  =  Z  (d2w/dx2)/[l+(dw/dx)2]3/2 

This  familiar  Green  Strain  beam  relationship  will  help  when  visualizing  what  is 
going  on  in  the  beams  however,  they  will  not  be  exactly  correct  as  the  program  of  choice, 
ABAQUS,  uses  a  different  method  to  calculate  its  strains.  ABAQUS  is  used  heavily  by 
BMW.  Therefore,  the  method  approach  ABAQUS  follows  is  very  useful  in  the  auto 
industry  which  deals  with  models  of  automobiles  crashing  where  rubber  and  polymers  are 
considered  as  materials.  ABAQUS  used  what  is  called  the  Green-Naghdi  (G-N)  rate. 

The  G-N  rate  can  make  use  of  the  usual  constitutive  relations  often  used,  but 
stress  is  related  to  a  rate  of  deformation  gradient  matrix.  That  is  what  makes  it  so  ideal 
when  we  have  large  movement.  In  addition,  since  it  is  time  dependent  it  can  be  used 
when  the  material  is  time  dependent  such  as  viscoelasticity,  viscoplasticity  as  well  as 
hyperelasticity.  It  thus  has  value  when  the  material  is  rubbery.  Therefore,  to  sum  it  up, 
ABAQUS  uses  a  large  displacement  time  dependent  theory  which  allows  rate  functional 
expressions  without  the  higher  order  terms.  It  removes  the  rigid  body  rotations  by  a 
corotational  system  of  coordinates  which  therefore  does  away  with  nonlinear 
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displacement  rate  gradients.  We  can  use  the  usual  constitutive  relationships  to  help 
visualize  what  is  happening,  but  keep  in  mind  that  the  stress  rate  is  related  to  the  rate 
deformation  tensor,  not  the  natural  strain. 

Because  ABAQUS/CAE/Explicit  follows  the  corotated  G-N  rate  relations,  a  load 
acting  vertically  down  at  the  end  of  a  cantilever  beam  will  translate  into  a  load  along  the 
axial  direction  and  in  the  normal  direction.  Thus,  as  the  rotation  gets  larger,  the  axial 
membrane  and  the  higher  order  nonlinear  membrane  terms  increase.  These  G-N  stress 
rate  terms  prove  that  the  analysis  developed  by  the  code  has  a  large  rotation  term  in  it. 

As  the  movement  gets  larger,  the  axial  loading  becomes  greater.  This  scenario  can  be 
expressed  very  nicely  with  the  Green  Strain  terms  as  shown  above. 

ABAQUS 

ABAQUS/CAE  will  be  used  for  the  all  Finite  Element  analysis  in  this  study.  The 
CAE  section  of  ABAQUS  will  be  used  to  handle  nonlinear  analysis  through  explicit 
integration. 

Geometric  Nonlinear  effects  in  ABAQUS  can  be  added  into  the  system  by  simply 
turning  the  NLGEOM  choice  ON.  Explicit  direct  integration  is  best  suited  to  wave 
propagation  problems  and  can  handle  nonlinearity  better  than  both  implicit  direct 
integration  and  the  modal  method.  Computer  storage  requirements  are  kept  to  a 
minimum  and  the  computational  cost  per  time  step  is  small.  The  drawback  is  the  critical 
time  step,  tcrit.  [3], 

An  explicit  step  in  ABAQUS  can  be  very  unstable.  If  At  or  tcrit  is  too  large, 
explicit  integration  fails;  if  At  is  unnecessarily  small,  the  time  it  takes  ABAQUS  to 
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compute  a  solution  can  take  a  very  long  time.  It  is  necessary  to  consider  both  when 
choosing  a  tcrit  [3].  ABAQUS  does  this  for  the  user  when  an  automatic  time  step  is 
selected. 

Frequency  Step 

There  are  many  types  of  natural  frequencies  of  a  beam  including  axial,  torsional, 
and  bending  directions.  This  study  will  concentrate  on  the  bending  natural  frequencies 
that  correlate  to  the  flapping  motion  of  a  beam.  The  bending  natural  frequencies  will  be 
in-plane  with  the  applied  force.  These  natural  frequencies  are  found  using  a  ‘Frequency 
Step’  in  ABAQUS/CAE.  The  natural  frequencies  are  matched  up  with  their  Eigenvectors 
to  ensure  they  are  the  correct  modes.  The  ABAQUS  results  should  match  the  theoretical 
bending  natural  frequency  equations.  Since  this  study  includes  both  fixed- free  and  fixed- 
fixed  beams,  both  are  shown  below. 

Theoretical  values  for  a  fixed-free  beam  bending  modes 

1.8752  I  El  ,  4.6902  [  El 

a),  = - J -  and  co-,  = - J -  (3) 

2 71  V  pAL4  In  V  z?AL4 

Theoretical  values  for  a  fixed-fixed  beam  bending  modes 
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Frequency  Response  Function 

A  Frequency  Response  Function  plots  a  magnitude  vs.  frequency  for  a  given  point 
in  a  system.  The  most  common  of  plotted  magnitudes  are  accelerations,  velocities,  and 
displacements.  An  example  of  a  common  damped  FRF  for  one  mode  is  shown  in  Figure 
4. 


16.00  20.00 

Frequency 


Figure  4.  Common  damped  FRF 

Because  this  Frequency  Response  Function  is  slightly  damped,  its  peak  Field  Data 
magnitude  does  not  approach  oo.  This  FRF  line  represents  a  single  location  on  the  beam 
where  data  was  acquired.  In  a  Finite  Element  simulation  this  location  will  be  on  any 
node  in  the  FEM  mesh.  In  laboratory  data,  these  locations  are  anywhere  that  there  is  data 
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taken  on  the  beam,  rather  it  be  taken  by  an  accelerometer,  laser  vibrometer,  or  high  speed 
photography. 

Frequency  Response  Functions  can  be  attained  in  various  ways  via  results  from 
experiments  in  a  laboratory  or  from  computer  simulations.  Frequency  sweeps  in  the 
frequency  domain  or  random  input  signals  in  the  time  domain  are  both  commonly  used. 

A  well  developed  FRF  allows  wn,  Wd,  wr  and  modal  damping  coefficients  t,„  to  be 
calculated,  while  multiple  FRFs  allow  Eigenvectors  to  be  computed. 

An  example  of  an  undamped  linear  Frequency  Response  Function  at  one  mode  is 
below  in  Figure  5.  This  FRF  represents  one  point  on  a  beam  subject  to  three  different 
increasing  loads,  FI,  F2  and  F2.  FO  represents  the  linear  natural  frequency  of  the  beam. 
In  a  linear  FRF  the  peak  for  any  given  mode  achieves  its  highest  point  at  the  value  of  the 
natural  frequency,  wn.  If  the  system  is  undamped,  the  response  will  approach  oo. 


Figure  5.  Undamped  Linear  FRF 
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An  example  of  an  undamped  nonlinear  Frequency  Response  Function  at  one 
mode  is  below  in  Figure  6.  This  FRF  once  again  represents  one  Finite  Element  (FE) 
node  on  a  beam  subject  to  two  different  loads  of  increasing  magnitude,  FI  and  F2.  FO 
represents  the  nonlinear  curve  that  curves  FI  and  F2  are  following.  This  curve  is  called 
the  Backbone  Curve  and  is  gone  over  in  more  depth  later.  In  a  nonlinear  FRF,  the  peak 
for  mode  n  can  attain  its  highest  point  at  or  near  con  if  the  structure  has  weak 
nonlinearity,  but  can  also  attain  its  highest  peak  at  an  co  either  less  or  greater  than  con  itself 
depending  on  the  structure’s  properties,  setup,  and  mode.  Like  the  undamped  linear  FRF, 
the  response  of  the  system  will  also  approach  co  if  the  system  is  truly  undamped.  If  a 
system  is  modeled  as  a  nonlinear  system,  its  magnitude  of  response  at  a  given  natural 
frequency  will  always  be  less  than  or  equal  to  the  same  system  modeled  linearly  [1], 


(Ref.  1) 

Figure  6.  Undamped  Nonlinear  FRF 
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To  produce  a  nonlinear  FRF  using  ABAQUS,  a  separate  run  is  accomplished  for 
each  individual  point  on  the  FRF  curve.  For  a  the  forcing  function  P  sin((ot),  P  will  be 
held  constant  while  co  is  varied.  Each  runs  steady-state  displacement  at  the  tip  is  recorded 
to  produce  a  particular  FRF. 

Rayleigh  Damping 

Damping  dissipates  energy  producing  the  envelope  in  which  a  free  vibration 
problem  will  decay  with  time.  It  also  limits  the  amplitude  of  a  Frequency  Response 
Function  and  prevents  response  from  going  to  oo.  One  way  to  estimate  damping  is  with 
Rayleigh  damping  or  proportional  damping.  In  Rayleigh  damping,  the  global  damping 
matrix,  [C]  is  proportional  to  both  the  mass  matrix  [M]  and  the  stiffness  matrix  [K]  by  the 
constants  on  and  a2.  [3], 

[C]  =  on  [M]  +  a2  [K]  (5) 

“The  ai[M]  contribution  damps  the  lowest  modes  most  heavily,  while  the  a2[K] 
contribution  damps  the  highest  modes  most  heavily.”  [3]  The  mass  dampening 
coefficient  a  is  found  by  the  equation 

oq  =  2-corco2-(^rco2-  £,2-coi)  /  (co22-  coi2)  (6) 

The  stiffness  dampening  coefficient  p  is  found  by  the  equation 

a2=  2-(%>-co2  -  §i-©i)  /  (co22  -  an2)  (7) 

where  can  and  are  the  nth  natural  in-plane  bending  frequencies  and  nth  modal  damping 
ratio  respectively.  For  all  experiments  in  this  study,  values  of  op  and  co2  are  found  using 
ABAQUS.  The  Finite  Elements  Models  are  built  in  ABAQUS  and  a  ‘Frequency  Step’ 
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set  up  to  calculate  the  ton  values.  The  c;n  values  are  simply  damping  ratios  that  can  be 
acquired  experimentally  using  the  half  power  method.  The  values  of  0.03,  0.003  and 
0.0003  were  all  assumed  in  the  previous  experiment  by  Ernst  and  Palazotto  and  the  E,n 
value  of  0.003  was  found  to  produce  the  most  accurate  results  to  the  Gordon,  Holkamp 
and  Spottswood  laboratory  experiments.  For  all  experiments  in  this  study,  c;n  is  assumed 
to  be  0.003. 

Free  Vibration  /  Log  Decrement  Step 

In  ABAQUS  a  free  vibration  of  an  undamped  linear  or  nonlinear  system  will 
oscillate  forever  without  losing  any  energy.  The  magnitude  of  the  Amplitude  Envelope  is 
the  same  at  time  t  =  100  as  it  is  at  t  =  1.  Once  damping  is  introduced  into  the  system,  the 
amplitude  of  the  oscillations  will  decrease  with  an  exponential  decay  following  the 
equation 

Amplitude  Envelope  =  Aenveiope  =  Ao*e^wt  =  A(t)  (8) 

The  log  decrement  step  will  produce  data  that  will  supply  both  the  Amplitude 
Envelope  and  the  Instantaneous  Frequency  for  development  of  the  backbone  curve. 

Backbone  Curve 

At  mode  n,  the  Backbone  Curve  represents  the  nonlinearity  of  the  structure  as  at  a 
given  point.  It  is  the  theoretical  line  that  the  FRF  follows.  Figure  7  shows  a  backbone 
curve  in  both  a  linear  and  nonlinear  system.  In  a  linear  system,  the  backbone  curve  is  a 
vertical  line  at  ron.  In  a  nonlinear  system,  the  backbone  curve  can  be  a  line  or  curve  that 
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starts  at  the  linear  value  of  con  and  leans  to  either  the  left,  softening,  or  to  the  right, 
hardening. 


Figure  7.  Backbone  Curve:  Linear  and  Nonlinear 


Figure  8  shows  a  damped  backbone  curve  experiencing  nonlinear  hardening.  The 
Frequency  Response  Function  has  a  definite  value  at  its  peak.  The  peak  is  at  a  higher 
value  of  co  than  the  linear  natural  frequency  as  the  FRF  follows  the  backbone  curve 
asymptotically. 
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Figure  8.  Damped  FRF:  Nonlinear  Hardening 

The  backbone  curve  lies  between  the  lower  and  upper  level  of  stability.  The 
dotted  line  in  Figure  9  represents  the  unstable  region  of  the  nonlinear  Frequency 
Response  Function.  This  region  cannot  be  simulated  in  a  laboratory  experiment.  For 
example,  if  a  frequency  sweep  is  done  on  the  beam  from  the  lower  frequencies,  the 
response  will  follow  the  upper  boundary  until  the  peak  is  reached.  The  response  will  then 
fall  off  to  a  stable  region  at  that  same  frequency.  The  unstable  region  along  with  the 
stable  region  under  the  upper  boundary  will  not  be  captured.  The  lower  boundary  stable 
region  can  be  captured  by  doing  a  frequency  sweep  from  the  higher  frequencies  to  lower 
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frequencies,  but  will  again  jump  up  to  the  upper  boundary,  missing  the  unstable  region 
altogether. 


(Ref.  1) 

Figure  9.  Backbone  Curve:  Unstable  Region 

The  backbone  curve  can  be  found  via  Schudt’s  method  [9]  which  includes  Sine- 
sweep  testing  and  jump-down  points.  However,  this  is  difficult,  and  not  only  requires 
experimental  capabilities,  but  also  requires  a  lot  of  time.  Therefore,  all  backbone  data 
found  in  this  study  will  be  found  using  the  Hilbert  transform  and  free  response  decrement 
data. 

Hilbert  Transform 

The  Hilbert  transform  named  after  David  Hilbert  has  been  used  in  signal 
processing  for  years.  It  permits  a  signal  in  the  time  domain  to  be  analyzed  in  the 
frequency  domain.  It  does  this  by  giving  the  original  signal  both  real  and  imaginary 
parts.  It  is  defined  by  the  integral  transform: 
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(9) 


H[y(t)]  =  y(t)  =  1/ji  J  ppr'dr 


.(t-r) 

where  y(t)  is  the  Hilbert  transform  of  the  original  signal  y(t). 

The  Hilbert  transform  shifts  the  phase  of  the  signal  by  -nil  while  still  keeping  the 
same  envelope,  amplitudes  and  other  spectral  components.  An  example  of  a  phase  shift 
to  a  signal  is  shown  below  in  Figure  10  [7], 


Figure  10.  Signal  and  Hilbert  transform  in  Envelope 


Notice  the  phase  shift  while  the  signals  both  lie  in  the  envelope  A(t).  The 
frequencies  of  the  two  signals  remain  the  same  as  well. 


Figure  11:  Amplitude  Envelope 


20 


The  peaks  of  the  signal  at  each  oscillation  can  be  plotted,  and  as  a  whole  produce 
the  envelope  of  the  signal,  A(t).  The  envelope  can  be  as  simple  as  a  straight  line  in  an 
undamped  system,  or  more  complicated  from  an  unstable  or  random  signal,  as  can  be 
seen  in  Figure  11.  In  a  damped  system,  A(t)  is  directly  related  to  the  decay  rate  t,  [7]. 

The  envelope  of  the  signal  is  very  important  in  the  Hilbert  transform  and  the 
breakdown  of  the  signal  attributes.  The  instantaneous  amplitudes  and  their  corresponding 
instantaneous  frequencies  are  completely  dependent  on  the  domain  of  the  Amplitude 
Envelope.  In  an  undamped  system  the  envelope  will  be  a  constant  horizontal  line  and  in 
turn  the  instantaneous  amplitude  and  frequency  are  constant  values  as  well.  For  this 
reason,  a  simulation  of  a  damped  beam  undergoing  free  vibration  is  perfect  for  use  with 
the  Hilbert  transform  to  find  the  backbone  curve. 

Using  the  Hilbert  transform,  what  is  known  as  an  Analytical  Signal,  Y,  can  be 
produced  [7], 

Y  =  y(t)  +  i*y(t)  (10) 

where  y(t)  is  referred  to  as  the  quadrature  or  conjugate  function  of  y(t).  This  is 
similar  to  the  Euler  formula  elz  =  cos(z)  +  i  •  sin(z)  where  sin(z)  is  the  Hilbert  transform 
ofcos(z)  [7]. 

y(t)  is  simplified  by  separating  it  into  its  two  parts,  y(t)  =  A(t)  c(t);  A(t)  for  the 
envelope  and  c(t)  for  the  oscillations  as  shown  in  Figure  12,  y(t)  then  becomes  y(t)  = 
A(t)c(t).  The  analytical  signal  becomes  Y(t)  =  A(t)[c(t)  +  i-  c  (t)]  [6]. 
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Figure  12:  Damped  Amplitude  Envelope 


By  breaking  up  Y(t)  into  smaller  sections  the  size  of  (t-to),  the  instantaneous 
amplitude,  instantaneous  phase,  and  instantaneous  frequencies  can  be  all  be  found. 

The  instantaneous  amplitude  can  be  found  by  the  equation  A(t)  =  |  Y(t)  | . 

The  instantaneous  phases  are  found  by  using  the  MATLAB  function  Phase  = 
angle[Y(t)]  [6], 

The  instantaneous  frequencies  can  be  found  by  taking  the  slope  of  the 
instantaneous  phase  over  the  same  time  period  (t-to)  [7]. 

This  transform  allows  somewhat  complex  systems  to  be  analyzed  in  the  time 
domain  and  allows  one  to  see  the  instantaneous  attributes  of  frequency,  phase,  and 
amplitude  [7].  These  values  will  be  used  in  the  estimation  of  the  Backbone  Curve  for  all 
displacements. 
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Duffing  Equation 

A  simple  model  used  to  approximate  the  nonlinearity  in  a  beam  is  the  SDOF 
Duffing  Oscillator 

y  +  C  y  +  coq^  y  +  Py^  =  P  sin(a)t)  (11) 

The  Duffing  equation  is  very  useful  as  it  simplifies  the  backbone  curve.  The 
backbone  curve  is  estimated  by 

Q2  =  co02  +  %  PA2  (12) 

Q  -  instantaneous  nonlinear  frequency 

A  -  instantaneous  displacement  amplitude 

co0  -  linear  natural  frequency  is  held  constant 

P  -  cubic  stiffness  coefficient 

In  the  equation  the  value  of  coo,  linear  natural  frequency  is  held  constant.  As  seen 
in  Figure  13,  this  coo  is  the  value  where  the  Hilbert  transform  data  curve  fit  crosses  the  Q 
axis  when  the  instantaneous  amplitude  is  0  [6].  The  Duffing  equation  is  meant  to  model 
a  1-D  nonlinear  spring-mass  system.  In  this  study  it  is  applied  to  the  beams  at  a  node 
which  undergoes  the  largest  displacement.  Each  node  on  the  beams  will  have  its  own 
unique  backbone  curve.  The  P  value  represents  the  strength  of  the  nonlinear  effect  on  the 
beam. 
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Figure  13:  Finding  Linear  Natural  Frequency 


MATLAB 

MATLAB  is  used  to  acquire  a  backbone  curve  by  utilizing  the  Hilbert  transform 
and  Duffing  equation  from  the  results  produced  by  ABAQUS/CAE  for  each  of  the  five 
Finite  Element  Models.  The  MATLAB  code  used  was  provided  by  Ernst  and  Palazotto 
[6], 

Research  Focus 

This  research  is  focused  on  the  evaluation  of  using  ABAQUS/CAE  to  solve 
nonlinear  problems  for  a  flapping  beam  and  whether  or  not  ABAQUS/CAE  results  used 
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with  MATLAB  via  the  Hilbert  transform  correctly  characterize  the  nonlinearity  of  beams 
based  on  previous  experiments.  Results  will  also  be  compared  to  one  another  to  find 
trends  of  nonlinearity  in  the  beams. 
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III.  Results  and  Discussion 


Chapter  Overview 

The  purpose  of  this  chapter  is  to  both  characterize  and  discuss  the  results  of  the 
beam  models  used  in  this  study.  The  subsequent  paragraphs  describe  the  five  beams. 

The  first  two  models  set  up  beams  with  clamped-clamped  supports  to  see  the  affects.  The 
last  three  models  set  up  beams  with  clamped-free  ends  to  depict  something  similar  to  a 
Micro  Air  Vehicle  wing.  All  of  the  beams  are  compared  in  their  effects  of  nonlinearity. 

Test  Subjects 

Clamped-Clamped  Steel  Beam  -  Beam  1 

To  start  out  the  study,  there  is  a  desire  to  develop  a  method  that  could  represent 
some  nonlinear  characteristics  in  a  relatively  timely  and  repeatable  process.  This  had 
been  done  before  by  Gordon,  Hollkamp  and  Spottswood  [9]  and  Ernst  [6]  for  the  same 
9.00”  clamped-clamped  Steel  beam.  The  study  accomplished  in  [9]  utilized  experimental 
data  with  the  Hilbert  transform  and  Duffing  equation.  The  study  in  [6]  used  Finite 
Element  analysis  data  along  with  the  Hilbert  transform  and  Duffing  equation.  Using  this 
same  beam  setup  will  both  verify  the  methods  used  in  this  study  as  well  as  give  a  baseline 
beam  model  to  compare  other  results  against.  This  beam  setup  will  be  referred  to  as 
Beam  1. 

Using  ABAQUS/CAE,  a  finite  element  model  was  developed  for  a  Clamped- 
Clamped  9”  long  Steel  beam  shown  in  Figure  15.  The  beam  itself  is  modeled 
symmetrically  as  a  4.5”  beam  in  two-dimensional  planar  modeling  space  in  ABAQUS  as 
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a  deformable  wire  in  the  1 -direction.  In  ABAQUS,  a  section  is  a  cross  sectional  area 
with  properties  and  dimensions  that  can  be  applied  to  wires.  A  section  is  applied  to  the 
wireframe  with  the  following  Steel  properties  and  dimensions: 


Figure  14.  Section  Dimensions:  Clamped-Clamped 


Table  1.  Section  Properties:  C-C  Steel  Beam 


t,  thickness  in  the  2-direction  =  0.031” 


w,  width  in  the  3-direction  =  0.5” 


E,  Young’s  Modulus  of  Elasticity  =  29.7  Mpsi 


n,  Poisson’s  ratio  =  0.33 


G,  Shear  Modulus  =11.1  Mpsi 

- 2 - 4~ 

p,  mass  density  =  7.36e-4  lbs  in" 
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A  cantilevered  or  “Encastre”  boundary  condition  is  applied  to  one  end  of  the  beam,  x  = 
0”.  A  Y-symmetry  boundary  condition  is  applied  at  x  =  4.5”.  A  load,  P,  is  also  applied 
to  the  midspan  of  the  beam  at  x  =  4.5”  as  shown  in  Figure  15. 


Figure  15.  C-C  Steel  beam 


Retrieving  the  Natural  Frequencies 

A  Frequency  Step  is  applied  in  ABAQUS  to  find  natural  bending  frequencies  of 
the  beam.  ABAQUS  shows  the  Eignevectors  along  with  the  natural  frequencies.  This 
makes  it  easy  to  separate  the  bending  natural  frequencies  from  the  torsional,  axial  and  out 
of  plane  bending  natural  frequencies. 
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Figure  16.  Frequency  step  picture 

The  results  of  the  frequency  step  are  shown  in  Table  2.  The  Eigenvectors  are  shown  in 
Figure  17  -  Figure  19. 


Figure  17.  First  Mode  Eigenvector;  wl  =  79.02  Hz 
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Figure  18.  Second  Mode  Eigenvector;  w2  =  217.70  Hz 


Figure  19.  Third  Mode  Eigenvector;  w3  =  426.91 


Table  2.  Natural  Frequencies:  Clamped-Clamped  Steel  Beam 


Theoretical  results 

00 1  =  79.12 

oo2  =  217.94 

ro3=  427.40 

ABAQUS  results 

00 1  =  79.02 

oo2  =  217.70 

ro3=  426.91 

Relative  Error 

0.13% 

0.11% 

0.11% 

The  Rayleigh  Damping  coefficients,  oq  and  oq,  are  solved  for  by  Equation  (6)  and 
(7)  and  an  estimated  t,i  =  '(^=  0.3%  along  with  the  ooi  and  (02  values  above.  The 
coefficients  oq  and  oq  are  placed  in  ABAQUS  in  the  advanced  damping  section  shown  in 
Figure  20. 

<xi  =  2.51 

—  note:  for  all  beams,  cq  «  oq  and  is  assumed  to  be  0 
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Figure  20.  Rayleigh  Damping  for  a  Section 

In  the  next  step  the  Rayleigh  damping  coefficients  will  allow  for  a  non  constant 
Amplitude  envelope.  The  damping  will  produce  a  log  decrement  decay  of  the  beams 
amplitude  over  time. 

Producing  the  Log  Decrement  Response 

Two  steps  are  necessary  to  produce  a  log  decrement  plot  in  ABAQUS.  Both  steps 
are  ‘Explicit,  General’  steps  with  the  nonlinearity  turned  on.  The  first  step  applies  a  load 
to  the  beam  to  cause  a  deflection  in  the  w  or  3-direction.  This  first  deflection  will  be  the 
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maximum  deflection  of  the  beam  throughout  both  time  steps.  In  the  first  step,  load  P  is 
set  as  a  ‘smooth  step’  function  with  an  maximum  value  of  0.5  lbf.  This  produces  a 
maximum  displacement  of  approximately  0.08  inches,  a  similar  deflection  to  past 
experiments  of  Ernst  and  Palazotto  [6].  The  load  is  applied  for  0.0125  seconds.  The 
smooth  step  allows  for  the  force  to  be  applied  with  a  zero  derivative  at  both  the  beginning 
and  end  of  the  0.0125  time  step.  A  generic  smooth  step  function  that  is  applied  to  this 
beam  as  well  as  the  others  is  shown  in  Figure  21. 


Figure  21.  Smooth  Step  Function 

In  the  second  step,  the  beam  is  allowed  to  freely  vibrate.  The  time  step  must  be 
long  enough  to  produce  an  adequate  amount  of  data  for  the  production  of  a  backbone 
curve.  The  load  from  step  1  is  deactivated.  Step  two  is  applied  for  length  of  time  of  1 .4 
seconds. 
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The  fixed-fixed  Steel  beam’s  free  vibration  response  is  below  in  Figure  22.  The 
beam  starts  out  at  a  displacement  of  0.07  inches  and  at  1.4  seconds  has  decayed  to  a 
displacement  of  0.015  inches.  During  this  decay  the  frequency  of  the  oscillation  is 
changing.  It  is  clearly  seen  that  the  oscillation  frequencies  at  the  beginning  of  the  decay 
are  higher  than  the  frequencies  towards  t  =  1.4  seconds.  This  plot  holds  all  of  the 
information  necessary  to  create  a  backbone  curve  of  the  first  mode  for  the  Fixed-Fixed 
Steel  beam.  This  data  is  saved  in  a  .rpt  file  and  opened  in  MATLAB. 
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Figure  22.  Free  Vibration:  Fixed-Fixed  Steel  Beam 

After  applying  the  Hilbert  transform  via  MATLAB,  the  free  vibration  signal  is 
broken  up  into  two  parts.  The  first  part  is  called  the  Amplitude  Envelope,  and  is  made  up 


Free  Vibration:  Fixed-Fixed  Steel  Beam 


33 


amplitude,  A,  (in) 


of  the  signals  maximum  amplitude  against  time.  The  Amplitude  Envelope  is  shown  in 
Figure  23.  The  second  part  is  the  called  the  Instantaneous  Frequency,  and  is  made  up  of 
the  signals  nonlinear  frequencies  at  any  given  time.  This  is  shown  in  Figure  24. 


Amplitude  Envelope:  Fixed-Fixed  Steel  Beam 


Figure  23.  Amplitude  Envelope:  Fixed-Fixed  Steel  Beam 
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Instantaneous  Nonlinear  Frequency:  Fixed-Fixed  Steel  Beam 


Figure  24.  Instantaneous  Frequency:  Fixed-Fixed  Steel  Beam 


These  two  parts  are  plotted  against  one  another  in  Figure  25  to  create  a  plot  of 
Amplitude  versus  Nonlinear  Frequency  in  what  essentially  is  the  part  of  the  backbone 
curve.  This  will  not  be  what  is  referred  to  as  the  Backbone  Curve  however  as  it  only 
contains  information  for  amplitudes  that  fell  within  the  range  of  the  free  vibration 
simulation  from  ABAQUS.  There  are  also  outliers  present  that  are  attributed  to  the 
beginning  and  ending  portions  of  the  ABAQUS  data. 
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amplitude,  A,  (in) 


Nonlinear  Frequency  vs  Amplitude:  Fixed-Fixed  Steel  Beam 


non-linear  freguency,  Omega,  (Hz) 

Figure  25.  Amplitude  vs  Nonlinear  Frequency:  Fixed-Fixed  Steel 

To  create  the  Backbone  Curve  for  all  amplitudes,  the  Amplitude  and  Nonlinear 
Frequency  values  are  squared  and  are  plotted  against  one  another  in  Figure  26.  The 
Estimated  Backbone  Curve  equation  (12)  is  applied  to  find  the  values  of  coo  and  p. 
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Amplitude2  \s  Omega2:  Fixed-Fixed  Steel  Beam 


Figure  26.  Duffing  Equation  Plot:  Fixed-Fixed  Steel  Beam 

The  following  information  was  obtained  for  the  first  mode  Backbone  Curve  for 
the  Fixed-Fixed  Steel  beam.  The  linear  frequency  ooo  =  79.02  Hz  and  the  coefficient  p  = 
1.814e8  in'2  s'2.  These  results  are  used  to  show  the  nonlinearity  of  the  beam  compared  to 
a  linear  backbone  curve  and  form  the  equation  Q2  =  C0()2  F  %  PA2. 

In  Figure  27  the  fixed-fixed  Steel  beams  backbone  curve  is  plotted  against  its  first 
natural  frequency  of  79.02  Hz.  The  backbone  curve  shows  that  the  beam  has  a  very 
strong  geometric  nonlinearity.  When  it  is  subject  to  a  steady  state  displacement  of  0.08 
inches,  the  nonlinear  natural  frequency  of  the  beam  more  than  doubles  from  79  Hz  to  165 
Hz.  This  change  in  the  Backbone  Curve  from  a  linear  to  nonlinear  analysis  is  clearly 
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seen  as  the  nonlinear  Backbone  Curve  is  showing  a  strong  hardening  affect  as  it  leans  to 
the  right.  This  can  be  attributed  to  the  increasing  contribution  that  the  membrane  strain 
term  has  in  the  nonlinear  Green  Strain  equation. 


"Steel  -  Fixed  Fixed  beam” 
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Figure  27.  Backbone  Curve:  Fixed-Fixed  Steel 


The  methods  used  in  this  study  show  close  results  with  previous  experiments  of 
Gordon,  Holkamp,  and  Spotswood,  and  Ernst  shown  in  Table  3 


Table  3.  Duffing  Equation  Results  Comparison 


Fixed-Fixed  Steel  Beams 

co0 

P 

GHS  experimental 

79.00 

1.80e8 

Ernst  ABAQUS 

79.0047 

1.84e8 

Current  ABAQUS 

79.02 

1.8140e8 
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Developing  the  Frequency  Response  Function  Plots 

As  mentioned  previously,  each  individual  point  on  the  Frequency  Response 
Function  is  found  by  applying  a  force  of  P  sin(wt)  at  the  midspan  of  the  beam.  ABAQUS 
then  simulates  the  beam  response  for  the  time  length  specified.  The  steady  state  response 
of  the  beam  is  then  measured  and  recorded.  For  the  Clamped-Clamped  Steel  beam  the 
following  individual  points  on  the  Frequency  Response  Functions  were  found. 


Force:  P  =  0.2lbf:  Frequency,  w  =  60-  150  Hz 

P  =  0.1:  w  =  60-135  Hz 
P=  0.025:  w  =  60-105 


The  Frequency  Response  Function  points  are  plotted  against  the  backbone  curve 
for  comparison.  Below,  Figure  28  shows  a  single  data  point’s  ABAQUS  run.  For  Beam 
1,  each  individual  data  run  takes  approximately  30  minutes  for  ABAQUS  to  simulate. 
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Figure  28.  Steady  State:  Fixed-Fixed  Steel 
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All  of  these  points  were  plotted  together  to  create  three  separate  Frequency 
Response  Curves.  These  curves  further  validate  the  Backbone  Curve  data  as  they  follow 
it  asymptotically  just  as  theory  predicts. 


"Steel  -  Fixed  Fixed  beam" 


Freq.  w.  hz 

Figure  29.  Frequency  Response  Function  Data:  Fixed-Fixed  Steel 


Clamped-Clamped  Aluminum  Beam  -  Beam  2 

Next  it  was  desired  to  look  at  a  more  flexible  beam  of  the  same  dimensions  and 
boundary  conditions.  The  same  beam  setup  as  Beam  1  is  used  while  only  changing  the 
material  properties  from  Steel  to  Aluminum.  The  main  values  that  changed  are  Young’s 
Modulus,  E  and  material  density,  p.  The  value  of  E  went  from  Esteei  =  29.7  Mpsi  to  Eaiuih 
=  10.6  Mpsi;  it  became  more  flexible  by  a  factor  of  almost  3.  The  density  of  the  beam 
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went  from  psteei  =  .000736  to  pAium  =  .000258;  it  became  less  dense  by  a  factor  of  almost 
3.  The  specific  Modulus  of  Elasticity,  E/p  only  changed  slightly  from  4.03el0  to 
4.09el0.  This  beam  will  show  how  a  more  flexible  Aluminum  beam  behaves  nonlinearly 
compared  to  a  stiffer  Steel  beam.  This  beam  setup  will  be  referred  to  as  Beam  2. 

The  clamped-clamped  model  for  the  Aluminum  beam  is  the  same  model  used  for 
Beam  1,  switching  appropriate  material  properties  shown  in  Table  4 
Table  4.  Section  Properties:  C-C  Aluminum  Beam 

t,  thickness  in  the  2-direction  =  0.031” 
w,  width  in  the  3-direction  =  0.5” 

E,  modulus  of  elasticity  =10.7  Mpsi 

u,  poisson’s  ratio  =  0.28 

G,  Shear  Modulus  =  4.14  Mpsi 
p,  mass  density  =  2.588e-4  lb-s2-in"4 


Retrieving  the  Natural  Frequencies 

A  Frequency  Step  is  applied  in  ABAQUS  to  find  natural  bending  frequencies  of 
the  beam.  The  natural  frequencies  of  the  Aluminum  beam  turned  out  to  be  just  slightly 
higher  than  the  Steel  beam  due  to  the  slightly  higher  specific  Modulus. 


Table  5.  Natural  Frequencies:  Fixed-Fixed  Aluminum  Beam 


Theoretical  results 

co  i  =  79.71 

ro2=  219.56 

(03=  430.58 

ABAQUS  results 

coi  =  79.60 

(02=  219.3 

(03=  429.9 

Relative  Error 

0.13% 

0.11% 

0.15% 
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The  Rayleigh  Damping  coefficients,  ai  and  0,2,  are  solved  for  using  the  newly 
found  00 1  and  002-  The  Rayleigh  Damping  coefficient  ai  =  2.202,  ct2  =  0.  The  coefficients 
ai  and  0(2  are  placed  into  ABAQUS  in  the  advanced  damping  section. 

Producing  the  Log  Decrement  Response 

Two  ‘Explicit,  General’  steps  are  setup  in  ABAQUS.  The  nonlinear  geometry 
option  is  turned  on.  In  the  first  step,  a  smooth  step  function  load  P  is  applied  in  the  w  or 
3-direction  with  a  maximum  value  of  0.5  lbf.  This  force  was  chosen  to  deflect  the  beam 
to  an  initial  displacement  of  0.07  inches,  the  same  amount  as  Beam  1  and  past 
experiments.  The  load  is  applied  for  0.0125  seconds.  In  the  second  step,  the  load  is 
deactivated  and  the  beam  is  allowed  to  freely  vibrate  until  the  deflections  were  about  0.01 
inches  at  their  maximum.  This  turned  out  to  be  for  1.6  seconds.  The  fixed-fixed 
Aluminum  beam’s  free  vibration  response  is  below  in  Figure  30.  The  beam  starts  out  at  a 
displacement  of  0.07  inches  and  at  1.6  seconds  has  decayed  to  a  displacement  of  0.01 
inches.  Just  like  the  fixed-fixed  Steel  beam,  the  changes  in  frequencies  are  visible  from 
the  start  to  finish  of  the  free  vibration  data.  The  data  saved  into  a  .rpt  file  and  is  opened 
in  MATLAB  for  backbone  curve  analysis. 
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displacement,  d,  (in) 


Free  Vibration:  Fixed-Fixed  Aluminum  Beam 


Figure  30.  Free  Vibration:  Fixed-Fixed  Aluminum  Beam 

MATLAB  and  incorporating  the  Hilbert  transform  yield  two  plots.  The 
Amplitude  Envelope  is  shown  in  Figure  31.  The  Instantaneous  Frequency  versus  time  is 
shown  in  Figure  32. 
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amplitude,  A,  (in) 


Amplitude  Envelope:  Fixed-Fixed  Aluminum  Beam 


Figure  31.  Amplitude  Envelope:  Fixed-Fixed  Aluminum  Beam 
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Instantaneous  Nonlinear  Freqnuency:  Fixed-Fixed  Aluminum  Beam 


Figure  32.  Instantaneous  Frequency:  Fixed-Fixed  Aluminum  Beam 


These  two  parts  are  plotted  against  one  another  as  in  the  previous  beam  in  Figure 
33  to  create  a  plot  of  Amplitude  versus  Nonlinear  Frequency.  Once  again,  there  are 
outliers  present  that  are  attributed  the  beginning  and  ending  portions  of  the  ABAQUS 
data. 
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amplitude,  A,  (in) 


Nonlinear  Frequency  vs  Amplitude:  Fixed-Fixed  Aluminum  Beam 


non-linear  freguency,  Omega,  (Hz) 

Figure  33.  Amplitude  vs  Nonlinear  Frequency:  Fixed-Fixed  Aluminum 

This  data  is  once  again  used  to  produce  a  straight  line  plot  of  Amplitude  and 
Nonlinear  Frequency  values  squared  as  shown  in  Figure  34.  The  Estimated  Backbone 
equation  (12)  is  applied  to  find  the  values  of  ooo  and  (3. 
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Amplitude2  vs  Omega2:  Fixed-Fixed  Aluminum  Beam 


amplitude2,  A2,  (in)2  x  10 3 

Figure  34.  Duffing  Equation  Plot:  Fixed-Fixed  Aluminum  Beam 

The  following  information  was  obtained  for  the  first  mode  Backbone  Curve  for 
the  Fixed-Fixed  Aluminum  beam.  The  linear  frequency  ooo  =  79.6  Hz  and  the  coefficient 
(3  =  1.84e8  in'2  s'2.  These  values  are  used  in  the  Estimated  Backbone  Curve,  Q2  =  coo2  +  % 
PA2.  In  Figure  35  the  fixed-fixed  Aluminum  beams  Backbone  Curve  is  plotted  against  its 
first  natural  frequency  of  79.6  Hz.  The  results  are  very  similar  to  the  fixed-fixed  Steel 
beam. 
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Figure  35.  Backbone  Curve:  Fixed-Fixed  Aluminum 


To  see  any  real  difference  between  the  Backbone  Curves,  they  must  be  plotted 
against  one  another  as  shown  in  Figure  36.  Both  beams  show  an  identical  shape 
however,  for  a  fixed-fixed  beam  with  similar  dimensions,  it  is  noticeable  that  a  more 
flexible  beam  exhibits  less  of  a  nonlinear  Backbone  Curve.  A  more  flexible  beam  will  be 
more  prone  to  bending  and  its  strains.  This  leads  to  a  higher  contribution  of  bending 
strain  in  the  total  strain,  a  lower  contribution  of  membrane  strains,  and  less  of  a  nonlinear 
result.  This  result  can  also  be  noticed  in  the  magnitude  of  the  (3  values  in  Table  12  where 
the  Aluminum  fixed- fixed  beam  has  a  lower  p  value  than  the  Steel  fixed- fixed  beam. 
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"Backbone  Curve  Comparison  -  Fixed  Fixed  beams: 
Steel  vs  Aluminum" 


Figure  36.  Backbone  Comparison:  Fixed-Fixed 


Clamped-  Aluminum  Beam  -  Beam  3 

This  is  the  first  beam  setup  that  will  have  a  clamped-free  boundary  condition 
resembling  the  setup  of  a  Micro  Air  Vehicle  wing.  Keeping  the  material  properties  the 
same  as  Beam  2,  nonlinear  characterization  can  be  compared  to  Beam  2  based  on  both 
beam  dimensions  and  beam  boundary  conditions.  This  beam  has  a  length  of  10.125”  and 
has  a  thickness  t  =  0.0625,  equal  to  two  times  the  thickness  of  Beams  1  or  2.  This  beam 
will  be  referred  to  as  Beam  3. 

Using  ABAQUS/CAE,  a  finite  element  model  was  developed  for  a  Clamped 
10.125”  long  Aluminum  beam.  Like  Beam  1  and  2,  this  beam  is  also  modeled  in  2 
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dimensional  planar  modeling  space  in  ABAQUS  as  a  10.125”  deformable  wire  in  the  1- 
direction.  A  section  is  applied  to  the  wireframe  with  the  following  Aluminum  properties 


and  dimensions: 


Figure  37.  Section  Dimensions:  Clamped  Aluminum 
Table  6.  Section  Properties:  C  Aluminum 

t,  thickness  in  the  2-direction  =  0.0625” 
w,  width  in  the  3-direction  =  0.5” 

E,  Young’s  Modulus  of  Elasticity  =  10.7  Mpsi 
n,  Poisson’s  ratio  =  0.28 
G,  Shear  Modulus  =  4.14  Mpsi 
p,  mass  density  =  2.588e-4  lb-s2-in'4 


A  cantilevered  boundary  condition  is  applied  to  one  end  of  the  beam,  x  =  0”  and  a 
load  P  is  applied  at  x  =  10.125”  as  shown  in  Figure  38 


50 


Figure  38.  C  Aluminum  Beam  setup 


Retrieving  the  Natural  Frequencies 

A  Frequency  Step  is  applied  to  find  the  natural  bending  frequencies  of  the  beam. 


Figure  39.  First  Mode  Eigenvector;  wl  =  19.911  Hz 


Figure  40.  Second  Mode  Eigenvector;  w2  =  124.49  Hz 
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Table  7.  Section  Properties:  Fixed  Aluminum  Beam 


Theoretical  results 

co  1  =  19.95 

ro2=  124.71 

ABAQUS  results 

coi=  19.91 

C02=  124.49 

Relative  Error 

0.2% 

0.17% 

The  Rayleigh  Damping  coefficients,  ai  and  0,2,  are  solved  for  using  these  new  coi 
and  C02  values.  The  Rayleigh  Damping  coefficient  ai  =  0.647  and  ct2  =  0.  The 
coefficients  ai  and  0.2  are  placed  ABAQUS  in  the  advanced  damping  section  to  allow  for 
a  damped  signal  to  be  simulated. 

Producing  the  Log  Decrement  Response 

Two  steps  are  made  to  produce  a  log  decrement  plot  in  ABAQUS.  As  before, 
both  steps  are  ‘Explicit,  General’  steps  with  the  nonlinear  geometry  option  turned  on. 

The  first  step  applies  a  load  to  the  beam  to  cause  a  deflection  in  the  w  or  3 -direction.  In 
the  first  step,  a  smooth  step  function  load  P  is  applied  with  a  maximum  value  of  0.7  lbf. 
This  leads  to  a  starting  deflection  of  4.5  inches.  The  same  smooth  step  type  of  load  from 
Beam  1  is  applied  for  0.0125  seconds. 

In  the  second  step,  the  beam  is  allowed  to  freely  vibrate.  The  load  is  deactivated 
and  the  beams  response  is  taken  for  12  seconds.  The  time  step  is  longer  than  the  previous 
two  examples  due  to  both  the  magnitude  of  the  initial  deflection  of  4.5  inches  compared 
to  0.08  inches  for  the  clamped-clamped  Steel  beam,  and  due  the  lower  damping 
coefficient  ai  estimated  by  a  lower  natural  frequency  of  this  beam  (19.9  Hz)  than  the 
previous  two  (80  Hz). 
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The  fixed-free  Aluminum  beam’s  free  vibration  response  is  below  in  Figure  41. 
The  beam  starts  out  at  a  displacement  of  4.5  inches  and  at  12  seconds  has  decayed  to  a 
displacement  of  0.5  inches.  The  log  decrement  step  was  taken  out  to  a  time  of  12 
seconds  in  order  to  allow  the  beam  oscillations  to  decay  down  to  an  amplitude  less  than  1 
inch.  If  the  beam  has  any  geometric  nonlinearity  present  then  there  will  be  some  change 
in  frequency  throughout  the  log  decrement  step.  This  change  in  frequency  is  not 
noticeable  as  it  was  with  the  fixed-fixed  beams,  however  the  Backbone  Curve  will 
provide  more  information.  The  data  saved  as  a  .rpt  file  and  opened  in  MATLAB  for 
backbone  curve  analysis. 


Free  Vibration:  Fixed-Free  Aluminum  Beam 


Figure  41.  Free  Vibration:  Fixed-Free  Aluminum  Beam 
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amplitude,  A,  (in) 


Using  MATLAB  and  incorporating  the  Hilbert  transform  yields  two  plots.  The 
Amplitude  Envelope  is  shown  in  Figure  42.  The  Instantaneous  Frequency  versus  time  is 
shown  in  Figure  43. 


Amplitude  Enwalope:  Fixed-Free  Alumimum  Beam 


Figure  42.  Amplitude  Envelope:  Fixed-Free  Aluminum 
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Instantaneous  Nonlinear  Frequency:  Fixed-Free  Aluminum  Beam 


Figure  43.  Instantaneous  Frequency:  Fixed-Free  Aluminum 


These  two  parts  are  plotted  against  one  another  in  Figure  44  to  create  a  plot  of 
Amplitude  versus  Nonlinear  Frequency.  Once  again,  the  outliers  present  are  attributed 
the  beginning  and  ending  portions  of  the  ABAQUS  data. 
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amplitude,  A,  (in) 


Nonlinear  Frequency  \s  Amplitude:  Fixed-Free  Aluminum  Beam 


non-linear  freguency,  Omega,  (Hz) 

Figure  44.  Amplitude  vs  Nonlinear  Frequency:  Fixed-Free  Aluminum 

The  Amplitude  (A)  and  Nonlinear  Frequency  (Q)  values  are  squared  and  are 
plotted  against  one  another  in  Figure  45.  The  estimated  Backbone  Curve  equation  (12)  is 
utilized  in  MATLAB  the  solve  for  the  values  of  coo  and  (3. 
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Amplitude2  vs  Omega2:  Fixed-Free  Aluminum  Beam 


Figure  45.  Duffing  Equation  Plot:  Fixed-Free  Aluminum 


The  following  information  was  obtained  for  the  first  mode  Estimated  Backbone 
Curve,  Q2  =  coo2  +  %  PA2,  for  the  Fixed-Free  Aluminum  beam.  The  linear  frequency  coo  = 
19.91 1  Hz  and  the  cubic  stiffness  coefficient  P  =  3.619  in'2  s'2.  This  is  much  lower  than 
the  Clamped-Clamped  p  values  of  1.8e8  in'2  s'2. 

In  Figure  46  Beam  3’s  Backbone  Curve  is  plotted  against  its  first  natural 
frequency  of  19.91 1  Hz.  The  beam  is  exhibiting  nonlinear  hardening  as  the  backbone 
curve  leans  to  the  right  of  the  linear  natural  frequency.  However,  the  backbone  curve  is 
not  showing  much  nonlinearity  at  all  just  like  the  frequency  of  oscillations  didn’t  seem  to 
be  changing  noticeably  in  the  log  decrement  plot  for  this  beam.  When  the  fixed- free 
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Aluminum  beam  experiences  steady-state  displacements  greater  than  10  times  the 
thickness  of  the  beam,  the  nonlinear  frequency  of  the  beam  changed  less  than  0.1  Hz. 


"Aluminum  -  Fixed  Free" 
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Figure  46.  Backbone  Curve:  Fixed-Free  Aluminum 


Due  to  the  fact  that  the  Backbone  Curves  for  the  clamped-clamped  Aluminum 
beam  and  the  clamped-free  Aluminum  beam  are  so  drastically  different,  if  they  are 
plotted  on  the  same  chart,  the  clamped-free  beam’s  Backbone  Curve  will  appear  to  be  a 
vertical  line.  The  fact  that  the  cubic  stiffness  coefficient  (P)  value  for  the  clamped- 
clamped  Aluminum  beam  is  about  7  orders  of  magnitude  larger  than  p  for  the  clamped- 
free  Aluminum  beam  shows  that  there  is  more  of  a  difference  in  the  beams  causing  the 
nonlinearity  than  either  their  thickness  or  length.  This  difference  can  only  be  attributed 
to  the  fact  that  one  beam  is  fixed  on  both  ends,  while  the  other  is  not.  As  shown  in  Figure 
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47  the  clamped-clamped  beam  has  two  walls  pulling  on  it  from  each  side.  As  the  force  P 
is  applied,  the  walls  constrain  the  beam  from  deflecting  or  bending  in  the  w  direction. 
This  is  not  only  causes  an  increase  in  the  membrane  strain  term  (du/dx)  but  also  a 
decrease  in  any  bending  strain  that  would  be  more  obvious  in  a  fixed-free  beam.  This 
adds  dramatically  to  the  value  of  the  cubic  stiffness  coefficient  P  for  a  fixed- fixed  beam. 


While  the  nonlinear  effects  on  the  clamped- free  Aluminum  beam  seems  small, 
especially  in  comparison  to  Beam  1,  other  beams  of  the  clamped- free  nature  will  be 
looked  at  to  bring  about  further  comparison. 

Clamped  Steel  Beam  -  Beam  4 

Next  it  was  desired  to  look  at  a  less  flexible  beam  of  the  same  dimensions  and 
boundary  conditions.  The  same  beam  setup  as  Beam  3  is  used  while  only  changing  the 
material  properties  from  Aluminum  to  Steel.  The  main  values  that  changed  are  Young’s 
Modulus,  E  and  material  density,  p.  The  value  of  E  went  from  EAiUm  =  10.6  Mpsi  to  Esteei 
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=  29.7  Mpsi;  it  became  less  flexible  by  a  factor  of  almost  3.  The  density  of  the  beam 
went  from  pAium  =  .000258  to  psteei  =  .000736;  it  became  more  dense  by  a  factor  of  almost 
3.  The  specific  Modulus  of  Elasticity,  E/p  only  changed  slightly  from  4.09el0  to 
4.03el0.  This  beam  will  again  show  how  nonlinear  characteristics  change  for  similar  a 
similar  Steel  and  Aluminum  beam  setup.  This  beam  setup  will  be  referred  to  as  Beam  4. 

The  clamped  Steel  beam  is  the  exact  model  used  for  Beam  3,  only  switching  the 
material  properties  found  below. 


Table  8.  Section  Properties:  C  Steel  Beam 


t,  thickness  in  the  2-direction  =  0.0625” 


w,  width  in  the  3-direction  =  0.5” 


E,  modulus  of  elasticity  =  29.7  Mpsi 


u,  poisson’s  ratio  =  0.33 


G,  Shear  Modulus  =11.1  Mpsi 
p,  mass  density  =  7.36e-4  lb-s2-in"4 


Retrieving  the  Natural  Frequencies 

A  Frequency  Step  is  applied  in  ABAQUS  to  find  natural  in-plane  bending 
frequencies  of  the  beam  found  in 


Table  9. 
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Table  9.  Section  Properties:  Fixed  Steel  Beam 


Theoretical  results 

00 1  =  19.80 

co2=  123.79 

ABAQUS  results 

(01=  19.76 

(o2=  123.56 

Relative  Error 

0.18% 

0.18% 

The  Rayleigh  Damping  coefficients,  ai  and  0,2  are  solved  for.  The  Rayleigh 
Damping  coefficient  ai  =  0.643  and  0.2  =  0.  The  coefficients  ai  and  ob  are  placed 
ABAQUS  in  the  advanced  damping  section. 

Producing  the  Log  Decrement  Response 

Two  ‘Explicit,  General’  steps  are  setup  in  ABAQUS.  The  nonlinear  geometry 
option  is  turned  on.  In  the  first  step,  a  smooth  step  function  load  P  is  applied  with  a 
maximum  value  of  1.0  lbf  that  deflects  the  tip  of  the  beam  to  4.5  inches.  The  load  is 
applied  for  0.0125  seconds.  In  the  second  step,  the  load  is  deactivated  and  the  beam  is 
allowed  to  freely  vibrate  for  12  seconds. 

The  fixed-free  Steel  beam’s  free  vibration  response  data  is  almost  identical  to  the  fixed- 
free  aluminum  beam  and  is  shown  below  in  Figure  48.  The  main  difference  is  the 
amount  of  force  it  takes  to  deflect  the  beam  to  the  same  initial  perturbation  of  4.5  inches. 
Any  frequency  change  is  not  seen  via  visual  inspection.  The  data  is  saved  as  a  .rpt  file 
and  opened  with  MATLAB  to  create  the  Backbone  Curve. 
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Free  Vibration:  Fixed-Free  Steel  Beam 


time,t,  (s) 


Figure  48.  Free  Vibration:  Fixed-Free  Steel  Beam 
The  Amplitude  Envelope  is  shown  in  Figure  49.  The  Instantaneous  Frequency 

versus  time  is  shown  in  Figure  50. 
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Amplitude  Envelope:  Fixed-Free  Steel  Beam 


Figure  49.  Amplitude  Envelope:  Fixed-Free  Steel 
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Instantaneous  Nonlinear  Frequency:  Fixed-Free  Steel  Beam 
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Figure  50.  Instantaneous  Frequency:  Fixed-Free  Steel 


These  two  parts  are  plotted  against  one  another  in  Figure  5 1  to  create  a  plot  of 
Amplitude  versus  Nonlinear  Frequency.  Once  again,  the  outliers  present  are  attributed 
the  beginning  and  ending  portions  of  the  ABAQUS  data. 
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amplitude,  A,  (in) 


Nonlinear  Frequency  vs  Amplitude:  Fixed-Free  Steel 


non-linear  freguency,  Omega,  (Hz) 

Figure  51.  Amplitude  vs  Nonlinear  Frequency:  Fixed-Free  Steel 

The  Amplitude  (A)  and  Nonlinear  Frequency  (Q)  values  are  squared  and  are 
plotted  against  one  another  in  Figure  52.  The  Estimated  Backbone  equation  Q2  =  coo2  +  % 
PA2  is  used  to  find  the  values  of  coo  and  p  by  fitting  it  to  the  data  points. 
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Figure  52.  Duffing  Equation  Plot:  Fixed-Free  Steel 


The  following  information  was  obtained  for  the  first  mode  Backbone  Curve  for 


the  Fixed-Free  Steel  beam.  The  linear  frequency  coo  =  19.73  Hz  and  the  coefficient  p  = 


3.951  in'2  s'2.  In  Figure  53  the  fixed-free  Steel  beams  backbone  curve  is  plotted  against 


its  first  natural  frequency  of  19.73  Hz. 


Because  the  p  value  for  Beam  4  is  very  similar  to  the  p  value  found  for  Beam  3, 


and  because  both  are  at  least  7  orders  of  magnitude  smaller  than  the  clamped-clamped 


beam’s  p  values,  it  is  apparent  that  clamped-clamped  beams  exhibit  much  greater 


nonlinear  tendencies  than  clamped-free  beams.  The  beam  setup  had  a  huge  influence  on 


the  nonlinear  characteristics  shown  in  both  the  Backbone  Curves  and  the  p  values.  The 


Amplitude* I 2 * * 5  vs  Omega2 
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boundary  condition  seems  to  be  the  key.  A  clamped-clamped  beam  has  much  more 
resistance  to  deflections  and  bending,  cutting  down  on  the  bending  strain  terms.  The 
resistance  is  caused  by  the  wall  boundary  conditions  pulling  on  the  beam  while  at  the 
same  time  causing  the  membrane  strain  terms  to  be  a  much  larger  contribution  of  the  total 
strain. 


"Steel  -  Fixed  Free  beam’ 
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Figure  53.  Backbone  Curve:  Fixed-Free  Steel 


Like  the  clamped-clamped  data,  the  results  for  the  fixed- free  Steel  beam  are 
similar  to  the  fixed- free  Aluminum  beam.  Figure  54  plots  the  two  backbone  curves 
against  one  another  for  the  two  fixed-fixed  beams.  Both  beams  are  completely  identical 
with  the  only  difference  being  the  difference  in  material  properties  between  Steel  and 
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Aluminum.  Similar  to  the  clamped-clamped  results,  this  plot  shows  that  for  a  fixed- fixed 
beam  with  similar  dimensions,  a  beam  made  of  Steel  has  a  slightly  stronger  nonlinear 
hardening  effect  than  a  more  flexible  Aluminum  beam.  This  result  can  also  be  noticed  in 
the  magnitude  of  the  cubic  stiffness  coefficient,  (3  values  in  Table  12.  The  Aluminum 
beams  both  have  lower  (3  values  than  their  corresponding  Steel  beams. 


"Backbone  Curve  Comparison  -  Fixed  Free  Beams: 


Non  Dimensionalized  Frequency.  w/wO.  (-) 


Figure  54.  Backbone  Comparison:  Fixed-Free:  Steel  vs  Aluminum 


Clamped  Composite  Beam  -  Beam  5 

The  previous  four  beams  are  all  isotropic,  and  while  their  information  is  useful  in 
characterizing  certain  nonlinear  effects,  Aluminum  and  Steel  may  not  be  likely 
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candidates  to  be  used  as  wing  material  for  a  future  Micro  Air  Vehicle  design.  A  material 
more  likely  to  be  used  would  be  some  sort  of  Composite  layup.  Therefore,  Beam  5  setup 
is  the  same  as  the  clamped-free  6.08  inch  Composite  beam  being  studied  by  Capt  Adam 
Tobias.  This  beam  is  the  thinnest  of  all  five  beams  at  t  =  0.02”. 

Using  ABAQUS/CAE,  a  finite  element  model  was  developed  for  the  Clamped 
6.08”  long  Composite  beam.  The  beam  is  a  4-ply  [0/90]s  layup  as  shown  in  Figure  55. 


Figure  55.  Composite  Beam  Ply  Layup 


The  composite  beam  is  modeled  in  3  dimensional  modeling  space  in  ABAQUS  as 
a  shell  with  a  total  thickness  of  .02”.  The  shell  is  divided  up  into  four  sections  or  lamina 
all  of  thickness  .005”.  The  length  of  the  beam  runs  parallel  with  the  1 -direction.  Each 
lamina  has  the  same  composite  properties  and  dimensions: 
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Table  10.  Section  Properties:  C  Composite  Beam 


tply,  thickness  of  each  ply  =  .005” 


w,  width  in  the  3-direction  =  0.5” 


El  =  19.2  Mpsi,  E2  =1.56Mpsi 


V12,  poisson’s  ratio  =  0.24 


G12  =  8.2e5,  G13  =  8.2e5,  G23  =  4.9e5 
p,  mass  density  =  .000163  lb-s2-in"4 


A  cantilevered  boundary  condition  is  applied  to  one  end  of  the  beam,  x  =  0”  and  a  load  P 
is  applied  in  the  2  or  3-direction  at  x  =  6.08”  as  shown  in  Figure  56. 


P 


Figure  56.  Composite  Beam  Load 

Retrieving  the  Natural  Frequencies 

A  Frequency  Step  is  applied  in  ABAQUS  to  find  the  in-plane  natural  bending 
frequencies  of  the  beam. 
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Figure  57.  First  Mode  Eigenvector;  wl  =  30.009  Hz 


Figure  58.  First  Mode  Eigenvector;  wl  =  188.48  Hz 


Table  11.  Section  Properties:  Fixed  Composite  Beam 


Theoretical  results 

001  =  30.00 

®2  =  188.00 

ABAQUS  results 

coi  =  30.009 

002=  188.48 

Relative  Error 

0.03% 

0.2% 

The  Rayleigh  Damping  coefficients,  ai  and  a.2,  are  solved  for  the  new  coi  and  002 
values.  The  Rayleigh  Damping  coefficient  ai  =  0.976  and  012  =  0.  The  ai  and  012  values 
are  placed  into  ABAQUS. 
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Producing  the  Log  Decrement  Response 

Two  steps  are  made  to  produce  a  log  decrement  plot  in  ABAQUS.  Both  steps  are 
‘Explicit,  General’  steps  with  the  nonlinear  geometry  option  turned  on.  The  first  step 
applies  a  load  to  the  beam  to  cause  a  deflection  in  the  3 -direction.  In  the  first  step,  a 
smooth  step  function  load  P  is  applied  with  a  maximum  value  of .  1  lbf.  The  load  is 
applied  for  0.0125  seconds  and  causes  a  maximum  deflection  of  3.7  inches  at  the  tip  of 
the  beam. 

In  the  second  step,  the  beam  is  allowed  to  freely  vibrate.  The  load  is  deactivated 
and  the  beams  response  is  taken.  The  load  from  step  1  is  deactivated.  Step  two  is  applied 
for  length  of  time  of  1.6  seconds.  The  fixed-free  Composite  beam’s  free  vibration 
response  is  below  in  Figure  59.  The  beam  starts  out  at  a  displacement  of  3.6  inches  and 
at  4  seconds  has  decayed  to  a  displacement  of  0.5  inches.  Like  the  other  fixed- free 
beams,  the  geometric  nonlinearity  of  the  beam  cannot  be  seen  visually  from  the  log 
decrement  plot.  The  Backbone  Curve  will  once  again  provide  more  information.  The 
data  is  saved  as  a  .rpt  file  and  opened  in  MATLAB  for  Backbone  Curve  analysis. 
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Free  Vibration:  Fixed-Free  [0/90]s  Composite 


Figure  59.  Free  Vibration:  Fixed-Free  Composite  Beam 
From  MATLAB  the  following  two  plots  were  determined.  The  Amplitude 

Envelope  is  shown  in  Figure  60.  The  instantaneous  frequency  versus  time  is  shown  in 
Figure  61. 
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amplitude,  A,  (in) 


Amplitude  Envelope:  Fixed-Free  [0/90]s  Composite 


Figure  60.  Amplitude  Envelope:  Fixed-Free  Composite 
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Instantaneous  Nonlinear  Frequency:  Fixed-Free  [0/90]s 


Figure  61.  Instantaneous  Freqeuncy  vs  Time:  Fixed-Free  Composite 


These  two  parts  are  plotted  against  one  another  in  Figure  62  to  create  a  plot  of 
amplitude  versus  nonlinear  frequency.  These  plots  look  very  similar  in  both  their  decay 
and  relatively  unchanging  frequencies  to  the  other  fixed  free  beams. 
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amplitude,  A,  (in) 


Nonlinear  Frequency  vs  Amplitude:  Fixed-Free  [0/90]s 
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Figure  62.  Amplitude  vs  Nonlinear  Frequency:  Fixed-Free  Composite 

The  Amplitude  (A)  and  Nonlinear  Frequency  (Q)  values  are  squared  and  are 
plotted  against  one  another  in  Figure  63.  The  values  of  coo  and  (3  are  determined  using  the 
Estimated  Backbone  equation  Q2  =  coq2  +  %  pA2. 
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Amplitude2  vs  Omega2 


Figure  63.  Duffing  Equation  Plot:  Fixed-Free  Composite 

The  following  information  was  obtained  for  the  first  mode  Backbone  Curve  for 
the  Fixed-Free  Composite  beam.  The  linear  frequency  coo  =  30.05  Hz  and  the  coefficient 
(3  =  47.45  in'2  s'2.  In  Figure  64  the  fixed-free  Composite  beams  Backbone  Curve  is 
plotted  against  its  first  natural  frequency  of  30.00  Hz.  The  plot  shows  that  just  like  all  of 
the  other  beams,  the  Composite  beam  is  undergoing  nonlinear  hardening. 


77 


Composite-  Fixed  Free " 


Freq.  w.  hz 

Figure  64.  Backbone  Curve:  Fixed-Free  Composite 


Figure  65  shows  the  fixed-free  Steel  and  Aluminum  Backbone  Curves  against  the 
fixed  free  Composite  curve.  All  backbone  curves  are  non-dimensionalized.  The 
displacements  are  divided  by  each  beams  length  while  the  nonlinear  frequencies  are 
divided  by  each  beams  linear  natural  frequency.  This  plot  shows  that  in  the  first  mode  a 
Composite  beam  experiences  a  higher  nonlinear  hardening  effect  than  either  an  isotropic 
Steel  or  Aluminum  beam  of  similar  dimensions.  This  shows  that  composite  beams  of 
[0/90]s  layups  are  more  affected  by  membrane  strain  than  the  isotropic  beams.  All  three 
clamped-free  beam  materials  show  very  weak  Backbone  Curves  compared  to  the 
clamped-clamped  beams.  When  the  clamped-free  beams  are  compared  to  one  another, 
the  isotropic  beams  Backbone  Curve  appears  even  weaker  than  it  was  before.  For 
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displacements  greater  than  20  times  the  thickness  for  isotropic  fixed- free  beams,  the 
nonlinear  characteristics  in  the  1st  mode  are  almost  negligible.  The  composite  beam  of 
similar  dimensions  shows  a  12  times  larger  contribution  of  the  membrane  strain  terms, 
‘/'-(dw/dx)  and  (du/dx).  This  result  can  also  be  noticed  in  the  magnitude  of  the  (3  values 
in  Table  12.  A  higher  positive  P  value  attributes  to  a  stronger  effect  of  nonlinear 
hardening. 

"Backbone  Curve  Comparison  -  Fixed  Free  Beams: 


Non  Dimensionalized  Frequency.  w/wQ.  (-) 

Figure  65.  Backbone  Comparison:  Fixed-Free:  Isotropic  vs  Composite 

It  is  not  just  the  composite  material  that  increases  the  nonlinearity  of  the  beam. 
The  same  composite  beam  dimensions  are  applied  to  a  composite  beam  with  a  [0/0]s 
layup.  The  Backbone  Curve  is  plotted  against  both  the  isotropic  Beam  3  and  4  Backbone 
Curves  as  well  as  the  Backbone  Curve  for  Beam  5’s  [0/90]s  layup.  This  is  shown  in 
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Figure  66  and  shows  that  the  actual  layup  of  the  Composite  beam  increases  the 
nonlinearity  effects  more  than  the  actual  values  of  Modulus  of  Elasticity  and  Density 
changing.  The  two  stronger  0  degree  layers  being  stacked  on  the  outside  of  the  beam 
allow  the  membrane  strain  to  be  more  of  a  contributing  factor  than  in  the  [0/0]s  beam. 


"Backbone  Curve  Comparison  -  Fixed  Free  Beams: 


Figure  66.  Backbone  Comparison:  Fixed-Free:  All 


Figure  67  compares  the  [0/90] s  Composite  with  a  density  =  0.000144  to  the  same 
[0/90]s  beam  with  a  slightly  higher  density  =  0.000163.  The  [0/0]s  beam  is  also 
compared.  This  shows  as  density  is  decreased,  the  membrane  strain  term  becomes  more 
of  a  contributing  factor,  and  the  nonlinear  effects  are  increased.  However,  it  is  not  as  big 
of  an  effect  as  changing  the  ply  layup. 


80 


"Backbone  Curve  Comparison  -  Fixed  Free  Beams: 


Non  Dimensionalized  Frequency.  w/wO.  (-) 

Figure  67.  Backbone  Comparison:  Fixed-Free  Composite  Beams 


We  already  saw  in  a  previous  comparison  that  the  clamped-clamped 
beams  experiences  a  much  stronger  effect  of  nonlinear  hardening  in  their  Backbone 
Curves  than  any  of  the  clamped-free  beams.  This  can  be  seen  in  both  the  p  values  both  in 
Table  12  and  visually  in  the  above  plots.  This  was  attributed  to  the  two  walls  in  the 
clamped-clamped  system  pulling  on  the  beam  and  adding  to  the  membrane  strain  effects. 
A  separate  Backbone  Curve  was  generated  for  three  different  nodes  on  the  [0/90]s 
Composite  beam  to  see  if  the  wall  pulling  on  the  beam  has  the  same  effect  in  a  clamped- 
free  system.  The  results  are  shown  in  Figure  68. 
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Figure  68.  Backbone  Comparison:  Fixed-Free  Composite:  3  Nodes 
Shown  in  Figure  69,  the  closer  the  point  on  the  beam  is  to  the  clamped  end  the 

less  of  an  effect  of  curviture  there  is.  This  corresponds  to  less  of  a  contribution  there  is  of 
the  bending  strain  term  (z  (d  w/dx  )).  This  increase  in  the  membrane  strain  contribution  . 
Likewise  the  point  on  the  free  end  of  the  beam  has  the  largest  curviture  and  the  largest 
contribution  of  bending  strain  terms.  For  this  reason,  the  cubic  stiffness  (3  value  is  always 
going  to  increase  as  a  point  on  a  fixed-free  beam  becomes  closer  to  the  fixed  end. 
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Figure  69.  Node  Placement  on  Composite  Beam:  1”,  3”  and  6.08”  from  fixed  end 


Chapter  Summary: 

The  five  beams  setups  show  how  both  boundary  conditions,  material  properties, 
and  beam  location  affect  the  nonlinear  characteristics  of  the  Backbone  Curve.  Table  12 

lists  all  of  the  Estimated  Backbone  Curve  values  for  Beams  1-5  from  the  equation  Q.  = 

2  2 

CGo  +  %  (3A  .  The  cubic  stiffness  coefficient  values  (3  are  a  direct  relationship  to  the 
strength  of  the  nonlinear  strength  of  each  beam. 
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Table  12.  Estimated  Backbone  Curve  Values 


Fixed-Fixed  Beams 

(Oo  -  Hz 

(3  -  in'2  s'2 

Steel  Beam  1 

79.09 

1.8140e8 

Aluminum  Beam  2 

80.07 

1.8105e8 

Fixed-Free  Beams 

N 

K 

1 

o 

s 

P  -  in  2  s'2 

Aluminum  Beam  3 

19.91 

3.619 

Steel  Beam  4 

19.73 

3.951 

Composite  [0/90] s  low  p 

30.00 

47.45 

Composite  [0/90]s  high  p 

28.28 

40.13 

Composite  [0/0]s  high  p 

30.01 

24.69 

In  fixed- fixed  beams,  the  vertical  load  P  is  quickly  translated  into  an  axial  force  as 
the  boundary  conditions  constrain  the  bending  of  the  beam.  This  leads  to  a  huge  p  value 
In  fixed-free  beams,  the  vertical  load  P  is  not  translated  into  an  axial  force  as 
quickly  as  in  a  fixed-fixed  beam.  It  takes  large  displacements  and  rotations  in  the  beam 
to  translate  the  P  force  into  an  force  in  the  axial  direction  of  the  beam.  However,  once 
this  begins  to  happen,  the  nonlinearity  of  the  fixed-free  beams  can  be  seen  in  the 
Backbone  Curves. 

The  less  flexible  Steel  beams,  Beam  1  and  Beam  4,  show  less  of  a  contribution  of 
bending  than  the  more  flexible  Aluminum  beams,  Beam  2  and  Beam  3.  This  leads  to 
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more  of  a  contribution  of  membrane  strain,  and  a  slightly  stronger  nonlinear  Backbone 
Curve  and  slightly  higher  p  values. 

The  [0/90]s  composite  beams  Backbone  Curve  leans  away  from  the  linear 
solution  even  more  so  as  the  P  value  is  12  times  higher  than  the  isotropic  beams.  The 
stronger  0  degree  plies  on  the  outide  of  the  beam  allow  for  a  higher  contribution  of 
membrane  strain,  and  thus  a  more  nonlinear  effect. 

Just  as  the  wall  boundary  condition  affected  the  fixed- fixed  beams  nonlinearity 
greatly,  as  a  point  on  a  fixed-free  beam  becomes  closer  to  its  fixed  end,  its  Backbone 
Curve  also  shows  a  much  higher  nonlinear  characteristic  with  much  higher  P  values. 
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IV.  Conclusions  and  Recommendations 


Chapter  Overview 

This  chapter  summarizes  the  important  findings  from  results  and  discussion.  The 
methods  used  in  finding  the  Backbone  Curve  using  ABAQUS/CAE  are  verified.  The 
trends  between  the  nonlinear  characteristics  and  the  following  are  reiterated: 

Beam  Setup 

Clamped-Clamped  vs  Clamped-Free 

Clamped-Free  beams  Backbone  Curve  for  different  points  on  beam 
-Each  point  has  a  new  Backbone  Curve 

Beam  Material 

Changes  in  Density,  p 
Changes  in  Stiffness,  E 
Changes  in  Composite  Beam  layup 
Changes  in  Thickness,  t 

Final  thoughts  on  how  this  affects  the  flapping  Micro  Air  Vehicle  study  are  made. 

Methods  Verified 

The  methods  used  in  this  study  showed  close  results  with  previous  experiments  of 
Gordon,  Holkamp,  and  Spotswood,  and  showed  a  slightly  better  result  than  the  previous 
ABAQUS  experiment  by  Ernst  as  shown  in  the  table  below 
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Table  13.  Past  and  present  Backbone  Curve  parameter  comparison 


Fixed-Fixed  Steel  Beams 

CDo 

P 

CHS  experimental 

79.00 

1.80e8 

Ernst  ABAQUS 

79.0047 

1.84e8 

Current  ABAQUS 

79.09 

1.8140e8 

Nonlinear  Characteristics  Summarized 

Clamped-Clamped  beams  exhibit  much  greater  nonlinear  tendencies  than 
Clamped-Free  beams.  The  beam  setup  had  a  huge  influence  on  the  nonlinear 
characteristics  shown  in  the  Backbone  curves.  The  cubic  stiffness  coefficient  values,  (3 
show  this  best  in  Table  12.  The  clamped-clamped  beams  both  have  (3  values  in  the  1.8e8 
range,  while  the  beta  values  for  the  clamped- free  beams  range  from  around  3  to  48.  The 
boundary  condition  seems  to  be  the  key.  A  clamped-clamped  beam  has  much  more 
resistance  to  deflections,  and  the  membrane  strain  term  is  a  bigger  player  because  of  the 
resisting  force  placed  on  the  beam  from  both  wall  boundary  conditions. 

The  closer  a  point  on  a  Clamped-Free  beam  is  to  the  wall  support,  the  greater  its 
nonlinear  tendencies.  This  is  shown  in  Figure  68  and  displays  that  membrane  strain 
contributions  are  highest  when  a  beam  has  a  wall  pulling  on  it.  The  Backbone  curve  for 
the  free  end  node  is  the  least  nonlinear  as  it  has  less  resistance  or  pull  on  it  axially  from 
the  wall. 
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A  more  flexible  beam  exhibits  less  of  a  nonlinear  Backbone  Curve.  Steel  beams, 
Beams  1  and  4,  show  slightly  higher  nonlinear  characteristics  effects  than  Aluminum 
beams,  Beams  2  and  3. 

As  thickness  of  a  beam  decreases,  the  membrane  strain  contribution  increases 
and  nonlinearity  increases.  Nonlinear  theory  mentions  that  the  thin  beams  experience 
more  contribution  through  their  membrane  than  thicker  beams.  No  studies  were  done  to 
verify  this,  however  the  composite  beam  was  3  times  thinner  than  the  isotropic  beams, 
and  experienced  a  higher  contribution  from  the  membrane  strain.  This  increase  in 
nonlinearity  was  due  to  both  the  thickness  of  the  beams  and  the  composite  material  layup. 

For  composite  beams,  a  [0/90] s  layup  is  more  affected  by  nonlinearity  than  a 
[0/0] s  layup.  This  is  shown  in  Figure  67.  It  shows  that  a  composite  beam  aligned  in 
somewhat  isotropic  manner  of  [0/0]  s  produces  a  similar  effect  as  an  isotropic  material  of 
Aluminum  or  Steel.  The  [0/90]  s  layup  has  produced  a  higher  contribution  of  membrane 
strain  leading  to  a  Backbone  Curve  that  falls  further  away  from  the  linear  solution. 

Composite  beams  are  more  affected  by  membrane  strain  and  the  nonlinear 
tendencies  that  come  with  it  than  beams  made  of  isotropic  materials.  This  is  shown  in 
Figure  65.  All  three  clamped-free  beam  materials  show  very  weak  Backbone  Curves 
compared  to  the  clamped-clamped  beams.  When  the  clamped-free  beams  are  compared 
to  one  another,  the  isotropic  beams  Backbone  Curve  looks  even  weaker.  There  doesn’t 
seem  to  be  nearly  as  much  membrane  strain  present  in  the  Aluminum  and  Steel  beams. 
That  being  said,  for  displacements  greater  than  20  times  the  thickness  for  isotropic  fixed- 
free  beams,  the  nonlinear  characteristics  in  the  1st  mode  are  negligible.  They  might 
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become  a  little  less  negligible  in  even  thinner  beams,  however  the  composite  beam  of  the 
same  thickness  would  still  show  a  larger  contribution  of  the  membrane  strain  terms, 
/^(dw/dx)2  and  (du/dx). 

Thoughts  on  Future  Research 

The  methods  prescribed  in  this  study  are  valuable  in  that  they  allow  us  to  find 
certain  nonlinear  characteristics  of  a  beam  setup  in  a  timely  manner  without  the  use  of 
laboratory  or  experimental  setup. 

While  the  clamped-clamped  beams  exhibit  extremely  high  nonlinear  effects,  they 
are  not  practical  models  when  representing  a  flapping  wing  Micro  Air  Vehicle.  When 
looking  at  the  first  mode  of  isotropic  clamped-free  beams,  their  nonlinearity  is  almost  if 
not  completely  negligible.  The  first  mode  nonlinearity  for  a  Composite  clamped-free 
beam  was  almost  7  times  as  noticeable.  Because  future  flapping  wing  MAVs  are  more 
likely  to  be  made  from  composite  materials,  the  nonlinearity  must  be  considered.  With 
wing  designs  that  will  certainly  have  a  much  smaller  thickness  and  density,  the 
membrane  stress  term  will  become  a  larger  contributor,  increasing  the  nonlinearity. 
Future  studies  will  be  needed  to  see  if  these  nonlinear  characteristics  will  be  a  factor. 
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Appendix 


Convergence:  Number  of  Elements  in  each  FEM 

The  accuracy  of  the  results  in  ABAQUS/CAE  relies  on  the  number  of  elements  in 
each  finite  element  model.  The  more  elements  there  are  in  the  model,  the  more  accurate 
the  analysis  will  be.  However,  with  the  increase  in  elements  comes  an  increase  in 
computational  expense.  For  each  model,  a  convergence  test  was  done  to  find  the  least 
number  of  elements  necessary  for  ABAQUS  to  find  wl  and  wl  within  .1%  of  their  actual 
values.  Each  models  true  natural  frequencies  were  found  using  at  least  n=1000  elements. 
The  1st  and  2nd  natural  bending  frequencies  were  then  found  for  the  models  using  n  =  1 
element,  increasing  n  until  the  values  converged  to  0.3%.  This  allows  the  ABAQUS 
model  to  keep  its  accuracy  while  cutting  down  on  the  computational  time  of  later  tests. 
The  Aspect  Ratio  of  each  element  is  also  taken  into  consideration.  Each  element  is 
preferred  to  have  an  Aspect  Ratio  of  1  to  1 .  The  number  of  elements  chosen  in  each 
Finite  Element  Model  was  chosen  to  also  keep  the  AR  close  to  a  1  to  1  ratio. 
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Clamped-Clamped  Steel  Finite  Element  Model 

ABAQUS  shows  that  the  Finite  Element  Model  only  needs  5  elements  to  keep 
accuracy  greater  than  0.3%  for  predicting  the  1st  and  3rd  natural  bending  frequencies. 

Five  elements  on  a  4.5”  beam  section  would  produce  5  elements  with  an  Aspect  Ratio  of 
almost  2  and  a  dimension  of  0.9”  x  0.5”.  Therefore  the  beam  was  seeded  in  ABAQUS  to 
have  9  elements  and  a  dimension  of  0.5”  x  0.5” 


Clamped-Clamped  Steel  Convergence:  Number  of  Elements  vs  1st  Natural  Frequency 


Figure  70.  Convergence  plots  C-C  Steel 
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Clamped-Clamped  Steel  Convergence:  Number  of  Elements  vs  3rd  Natural  Frequency 


n,  number  of  elements 

Figure  71.  Convergence  plots  C-C  Steel 
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Clamped  Aluminum  Finite  Element  Model 


ABAQUS  shows  that  the  Finite  Element  Model  only  needs  in  between  20 
elements  to  keep  accuracy  greater  than  0.3%  for  predicting  the  1st  and  2nd  natural  bending 
frequencies.  Twenty  elements  on  a  10.125”  beam  section  would  produce  20  elements 
with  an  Aspect  Ratio  of  1.0125  and  a  dimension  of  0.50625”  x  0.5”.  Therefore  the  beam 
was  seeded  in  ABAQUS  to  have  20  elements. 


Clamped  Aluminum  Convergence:  Number  of  Elements  vs  1st  Natural  Frequency 


Figure  72.  Convergence  plots  C-C  Aluminum 
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r,  2nd  natural  frequency,  Hz 


Clamped  Aluminum  Convergence:  Number  of  Elements  vs  2nd  Natural  Frequency 


Figure  73.  Convergence  plots  C  Aluminum 
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Clamped  Composite  Finite  Element  Model 

ABAQUS  shows  that  the  Finite  Element  Model  only  needs  around  80  elements  to  keep 
accuracy  greater  than  0.3%  for  predicting  the  1st  and  2nd  natural  bending  frequencies. 
Ninety  elements  on  a  0.5”  x  6”  beam  would  produce  90  elements  in  a  3  x  30  grid  with  an 
Aspect  Ratio  of  1 .2  and  a  dimension  of  0.2”  x  0. 166”.  The  beam  was  seeded  in 
ABAQUS  to  have  90  elements. 


Clamped  Composite  Convergence:  Number  of  Elements  vs  1st  Natural  Frequency 


Figure  74.  Convergence  plots  C  Composite 
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Clamped  Composite  Convergence:  Number  of  Elements  vs  2nd  Natural  Frequency 


Figure  75.  Convergence  plots  C  Composite 


96 


Clamped-Clamped  Aluminum  Finite  Element  Model 

Like  the  Clamped-Clamped  Steel  beam,  ABAQUS  shows  that  the  Finite  Element 
Model  only  needs  5  elements  to  keep  accuracy  greater  than  0.3%  for  predicting  the  1st 
and  3rd  natural  bending  frequencies.  Five  elements  on  a  4.5”  beam  section  would 
produce  5  elements  with  an  Aspect  Ratio  of  almost  2  and  a  dimension  of  0.9”  x  0.5”. 
Therefore  the  beam  was  seeded  in  ABAQUS  to  have  9  elements  and  a  dimension  of  0.5” 
x  0.5” 


Clamped-Clamped  Aluminum  Convergence:  Number  of  Elements  vs  1st  Natural  Frequency 


Figure  76.  Convergence  plots  C  C  Aluminum 
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Clamped-Clamped  Aluminum  Convergence:  Number  of  Elements  vs  3rd  Natural  Frequency 


n,  number  of  elements 

Figure  77.  Convergence  plots  C-C  Aluminum 
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Clamped  Steel  Finite  Element  Model 

Like  the  Clamped  Aluminum  beam,  ABAQUS  shows  that  the  Finite  Element  Model  only 
needs  in  between  20  elements  to  keep  accuracy  greater  than  0.3%  for  predicting  the  1st 
and  2nd  natural  bending  frequencies.  Twenty  elements  on  a  10.125”  beam  section  would 
produce  20  elements  with  an  Aspect  Ratio  of  1.0125  and  a  dimension  of  0.50625”  x  0.5”. 
Therefore  the  beam  was  seeded  in  ABAQUS  to  have  20  elements. 


Clamped  Steel  Convergence:  Number  of  Elements  vs  1st  Natural  Frequency 
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,  2nd  natural  frequency,  Hz 


Clamped  Steel  Convergence:  Number  of  Elements  vs  2nd  Natural  Frequency 


Figure  79.  Convergence  plots  C  Steel 
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MATLAB  CODE  for  Backbone  Curve  Estimation 


%BackboneCurveEstimation 
%Robert  Walker 

%%original  code  provided  by  Dr  Anthony  Palazottl  and  Matthew  Ersnt. 

clear  all; 

clc; 

%Block  1  DATA 


%Xel:  Clamped-Clamped  Steel  beam:  (.5  x  .031  x  9)  1.4  sec,  w  =  79.0  Hz 

%EalumX:  Clmpd-Clmpd  Alum  beam:  (.5  x  .031  x  9)  1.6  sec,  w  =  79.73  Hz 


%X12BA2 :  Clmpd-Free  Alum  beam  (.5  x  .0625  x  10.125)  12  sec,  w  =  19.91  Hz 
%CSteelX:  Clmpd-Free  Steel  beam  (.5  x  .0625  x  10.125)  12  sec,  w  =  19.76  Hz 

%XC442:  Clmp-Free  Comp  [ 0 / 0 ] s  beam  (.5  x  .02  x  6.08)  4  sec,  w  =  30.0  Hz 

%C0 90X1  CO 90X2  and  C090X3:  [0/90]s  den  =  0.166,  4  sec,  w  =  28  Hz 

%Cden090Xl  Cden090X2  and  Cden090X3:  [0/90] s  with  den  =  0.1445,  4  sec,  w  =  30  Hz 


load  '  (place. rpt  file  here)  ’ 

X  =  (place. rpt  file  here) 

omegal  =  30;  %place  natural  frequency  of  beam  here 
%Block  2 


[L  W]  =  size (X) ; 

interval  =  round ( 2 *L/X (L, 1 ) /omegal ) ;  %sets  up  individual  increments 
%interval  =  51 


%Block  3 

skip  =2;  %  sets  up  the  overlaping  of  the  data  in  each  increment 


%Block  4 

time  =  X (L, 1 )  -  X (1,1); 


%Block  5 

Z  =  hilbert (X ( : , 2 ) ) ; 

%Block  6 
i  =  0 

for  h=l : skip :  (L-2 ^interval); 


i  =  i  + 1  ; 
rotation  =  0; 
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%Block  7 

for  j =1 : interval ; 

Za ( j )  =  Z (h+j ) ; 

%Block  8 

absZa(j)  =  abs(Za(j)); 

%Block  9 

Zangleclean ( j )  =  angle (Za (j  ))  ; 

%Block  10 

Zangle(j)  =  angle (Za (j )) ; 
if  j==i; 

Zangle ( j ) ; 

elseif  Zangle (j)<=0  &  Zangleclean (j -1 ) >0 
rotation  =  rotation+2 *pi ; 

Zangle (j)  =  Zangle (j)  +rotation; 

else 

Zangle (j)  =  Zangle (j)+  rotation; 

end 

end 

b  =  ( 0 : 1 : interval-1 )  ./ (L/time) ; 

%Block  11 

PI  =  polyf it (b, log (absZa) , 1 )  ; 

%if  h==5*skip+l 
%f igure ( 6 ) 

%axis  tight 

%plot (b, log (absZa)  ,  ' *-  1  ) 

%end 

%Block  12 

A (i)  =  exp (PI  (2) )  ; 

%Block  13 

P2  =  polyf it (b, Zangle, 1 ) ; 

%if  h==5*skip+l 
%f igure ( 7 ) 

%axis  tight 
%plot (b, Zangle) 

%end 

%Block  14 

omega (i)  =  P2(l)/2/pi; 

end 

tine  =  time/L^interval/2+time/L^skip^ ( 1 : i ) ; 

%Block  15 
figure ( 1 ) 

plot (X ( : , 1) ,X  (  : , 2) ) 

xlabel ( ' time, t ,  (s) ') 

ylabel (’ displacement ,  d,  (in)1) 

%Block  16 
figure (2 ) 

%subplot (2,2,1) 

plot  (tine, omega,  ' r .  ’  ) 

xlabel ( ' time,  t,  (s) ') 

ylabel ( 1  non-linear  freqnecy.  Omega,  (Hz)  ') 
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%Block  17 
figure  ( 3 ) 

%subplot  (2,2,2) 
plot (tine, A, '  .  ’  ) 
xlabel ( ' time,  t,  (s) ') 
ylabel (' amplitude.  A,  (in)') 

%Block  18 
figure  ( 4 ) 

%subplot (2,2,3) 

plot (omega, A, ’ . ' ) 

ylabel (' amplitude,  A,  (in) ’) 

xlabel (' non-linear  freguency.  Omega,  (Hz) '  ) 

%Block  19 
figure  ( 5 ) 

%subplot  (2,2,4) 

plot (A. A2 , omega . A2 , 1 . ' ) 

xlabel (’ amplitudeA2 ,  AA2,  (in)A2') 

ylabel (’ non-linear  frequencyA2,  0megaA2,  (Hz)A2’) 
%Block  20 

P3  =  polyf it (A . A2 , omega . A2 , 1 )  ; 

%Block  21 

wO  =  sqrt (P3  (2 ) ) 

%Block  22 

Beta  =  4/3^P3 (1) * (2*pi) A2 
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